Changeset 2795


Ignore:
Timestamp:
01/12/10 09:36:38 (15 years ago)
Author:
Mathieu Morlighem
Message:

removed useless code

Location:
issm/trunk/src/c/Bamgx
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Bamgx/Mesh2.h

    r2794 r2795  
    889889  Vertex * NearestVertex(Icoor1 i,Icoor1 j) ;
    890890  Triangle * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
    891   void Write(const char * filename,const TypeFileMesh type = AutoMesh);
    892   void Write_am_fmt(std::ostream &) const ;
    893   void Write_am(std::ostream &) const ;
    894   void Write_ftq(std::ostream &) const ;
    895   void Write_nopo(std::ostream &) const ;
    896   void Write_msh(std::ostream &) const ;
    897   void Write_amdba(std::ostream &) const ;
    898 
    899   void Read(MeshIstream &,int version,Real8 cutoffradian);
     891
    900892  void ReadFromMatlabMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
    901893  void WriteBamgMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
    902   void Read_am_fmt(MeshIstream &);
    903   void Read_amdba(MeshIstream &);
    904   void Read_am(MeshIstream &);
    905   void Read_nopo(MeshIstream &);
    906   void Read_ftq(MeshIstream &);
    907   void Read_msh(MeshIstream &);
    908894
    909895  void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);

  • issm/trunk/src/c/Bamgx/MeshRead.cpp

    r2794 r2795  
    268268}
    269269
    270 void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian) {
    271 
    272 
    273         long int verbosity=0;
    274 
    275   Real8 hmin = HUGE_VAL;// the infinie value
    276   //  Real8 MaximalAngleOfCorner = 10*Pi/180;//
    277   Int4 i;
    278   Int4 dim=0;
    279   Int4 hvertices =0;
    280   Int4 ifgeom=0;
    281   Metric M1(1);
    282   if (verbosity>1)
    283     cout << "  -- ReadMesh " << f_in.CurrentFile  << " Version = " << Version << endl;
    284   int field=0;
    285   int showfield=0;
    286   while (f_in.cm())
    287     { 
    288       field=0;
    289       char  fieldname[256];
    290       if(f_in.eof()) break;
    291       f_in.cm() >> fieldname ;;
    292       if(f_in.eof()) break;
    293       f_in.err() ;
    294       if(verbosity>9)
    295         cout <<  "    " << fieldname << endl;
    296       if (!strcmp(fieldname,"MeshVersionFormatted") )
    297          f_in >>   Version ;
    298       else if(!strcmp(fieldname,"End"))
    299          break;
    300       else if (!strcmp(fieldname,"Dimension"))
    301          {
    302            f_in >>   dim ;
    303            assert(dim ==2);
    304          }
    305       else if  (!strcmp(fieldname,"Geometry"))
    306         {
    307           char * fgeom ;
    308           f_in >> fgeom;
    309           //      if (cutoffradian>=0) => bug if change edit the file geometry
    310           //  Gh.MaximalAngleOfCorner = cutoffradian;
    311           if (strlen(fgeom))
    312               Gh.ReadGeometry(fgeom);
    313           else
    314             {
    315               // include geometry
    316               f_in.cm();
    317               Gh.ReadGeometry(f_in,fgeom);
    318             }
    319 
    320           Gh.AfterRead();
    321           ifgeom=1;
    322           delete [] fgeom;
    323         }
    324       else if (!strcmp(fieldname,"Identifier"))
    325         {
    326           if (identity) delete [] identity;
    327           f_in >> identity;
    328         }
    329       else if (!strcmp(fieldname,"hVertices"))
    330          { hvertices =1;
    331            Real4 h;
    332            for (i=0;i< nbv;i++) {
    333             f_in >>  h ;
    334             vertices[i].m = Metric(h);}
    335          }
    336       else if (!strcmp(fieldname,"Vertices"))
    337          {
    338            assert(dim ==2);
    339            f_in >>   nbv ;
    340            if(verbosity>3)
    341              cout << "   Nb of Vertices = " << nbv << endl;
    342            nbvx=nbv;
    343            vertices=new Vertex[nbvx];
    344            assert(vertices);
    345            ordre=new (Vertex* [nbvx]);
    346            assert(ordre);
    347            
    348            nbiv = nbv;
    349            for (i=0;i<nbv;i++) {
    350              f_in >>
    351                          vertices[i].r.x >> 
    352                    vertices[i].r.y >>
    353              vertices[i].ReferenceNumber  ;
    354              vertices[i].DirOfSearch =NoDirOfSearch;
    355                    vertices[i].m=M1;
    356              vertices[i].color =0;}
    357            nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    358            triangles =new Triangle[nbtx];
    359            assert(triangles);
    360            nbt =0;
    361          }
    362       else if (!strcmp(fieldname,"Triangles"))
    363         {
    364           if (dim !=2)
    365             cerr << "ReadMesh: Dimension <> 2" <<endl , f_in.ShowIoErr(0);
    366           if(! vertices || !triangles  || !nbv )
    367             cerr << "ReadMesh:Triangles before Vertices" <<endl,
    368               f_in.ShowIoErr(0);
    369           int NbOfTria;
    370           f_in >>  NbOfTria ;
    371           if(verbosity>3)
    372             cout << "   NbOfTria = " << NbOfTria << endl;
    373           if (nbt+NbOfTria >= nbtx)
    374             cerr << "ReadMesh: We must have 2*NbOfQuad + NbOfTria  = "
    375                  << nbt + NbOfTria<<"  < 2*nbv-2 ="  << nbtx << endl,
    376               f_in.ShowIoErr(0);
    377           //      begintria = nbt;
    378           for (i=0;i<NbOfTria;i++) 
    379             {
    380               Int4 i1,i2,i3,iref;
    381               Triangle & t = triangles[nbt++];
    382               f_in >>  i1 >>  i2 >>   i3 >>   iref ;
    383               t = Triangle(this,i1-1,i2-1,i3-1);
    384               t.color=iref;
    385             }
    386           // endtria=nbt;
    387         }
    388       else if (!strcmp(fieldname,"Quadrilaterals"))
    389         {
    390 
    391           if (dim !=2)
    392             cerr << "ReadMesh: Dimension <> 2" <<endl , f_in.ShowIoErr(0);
    393           if(! vertices || !triangles  || !nbv )
    394             cerr << "ReadMesh:Quadrilaterals before Vertices" <<endl,
    395               f_in.ShowIoErr(0);
    396           f_in >>   NbOfQuad ;
    397           if(verbosity>3)
    398             cout << "   NbOfQuad= " << NbOfQuad << endl;
    399           if (nbt+2*NbOfQuad >= nbtx)
    400             cerr << "ReadMesh: We must have 2*NbOfQuad + NbOfTria  = "
    401                  << nbt + 2*NbOfQuad <<"  < 2*nbv-2 ="  << nbtx << endl,
    402               f_in.ShowIoErr(0);
    403           //  beginquad=nbt;
    404           for (i=0;i<NbOfQuad;i++) 
    405             {
    406               Int4 i1,i2,i3,i4,iref;
    407               Triangle & t1 = triangles[nbt++];
    408               Triangle & t2 = triangles[nbt++];
    409               f_in >>  i1 >>  i2 >>   i3 >> i4 >>    iref ;
    410               t1 = Triangle(this,i1-1,i2-1,i3-1);
    411               t1.color=iref;
    412               t2 = Triangle(this,i3-1,i4-1,i1-1);
    413               t2.color=iref;
    414               t1.SetHidden(OppositeEdge[1]); // two time  because the adj was not created
    415               t2.SetHidden(OppositeEdge[1]);
    416             }
    417           // endquad=nbt;
    418         }
    419       else if (!strcmp(fieldname,"VertexOnGeometricVertex"))
    420         {
    421            f_in  >> NbVerticesOnGeomVertex ;
    422            if(verbosity>5)
    423              cout << "   NbVerticesOnGeomVertex = "   << NbVerticesOnGeomVertex << endl
    424                   << " Gh.vertices " << Gh.vertices << endl;
    425            if( NbVerticesOnGeomVertex)
    426              {
    427                VerticesOnGeomVertex= new  VertexOnGeom[NbVerticesOnGeomVertex] ;
    428                if(verbosity>7)
    429                  cout << "   VerticesOnGeomVertex = " << VerticesOnGeomVertex << endl
    430                       << "   Gh.vertices " << Gh.vertices << endl;
    431                assert(VerticesOnGeomVertex);
    432                for (Int4 i0=0;i0<NbVerticesOnGeomVertex;i0++) 
    433                  {
    434                    Int4  i1,i2;
    435                    //VertexOnGeom & v =VerticesOnGeomVertex[i0];
    436                    f_in >>  i1 >> i2 ;
    437                    VerticesOnGeomVertex[i0]=VertexOnGeom(vertices[i1-1],Gh.vertices[i2-1]);
    438                  }
    439              }
    440          }
    441       else if (!strcmp(fieldname,"VertexOnGeometricEdge"))
    442          {
    443            f_in >>  NbVerticesOnGeomEdge ;
    444            if(verbosity>3)
    445              cout << "   NbVerticesOnGeomEdge = " << NbVerticesOnGeomEdge << endl;
    446            if(NbVerticesOnGeomEdge)
    447              {
    448                VerticesOnGeomEdge= new  VertexOnGeom[NbVerticesOnGeomEdge] ;
    449                assert(VerticesOnGeomEdge);
    450                for (Int4 i0=0;i0<NbVerticesOnGeomEdge;i0++) 
    451                  {
    452                    Int4  i1,i2;
    453                    Real8 s;
    454                    //VertexOnGeom & v =VerticesOnGeomVertex[i0];
    455                    f_in >>  i1 >> i2  >> s;
    456                    VerticesOnGeomEdge[i0]=VertexOnGeom(vertices[i1-1],Gh.edges[i2-1],s);
    457                  }
    458              }
    459          }
    460       else if (!strcmp(fieldname,"Edges"))
    461         {
    462           Int4 i,j, i1,i2;
    463           f_in >>  nbe ;
    464           edges = new Edge[nbe];
    465           if(verbosity>5)
    466             cout << "     Record Edges: Nb of Edge " << nbe << " edges " <<  edges << endl;
    467           assert(edges);
    468           Real4 *len =0;
    469           if (!hvertices)
    470             {
    471               len = new Real4[nbv];
    472               for(i=0;i<nbv;i++)
    473                 len[i]=0;
    474             }
    475           for (i=0;i<nbe;i++)
    476             {
    477               f_in  >> i1  >> i2  >> edges[i].ref ;
    478                  
    479               assert(i1 >0 && i2 >0);
    480               assert(i1 <= nbv && i2 <= nbv);
    481               i1--;
    482               i2--;
    483               edges[i].v[0]= vertices +i1;
    484               edges[i].v[1]= vertices +i2;
    485               edges[i].adj[0]=0;
    486               edges[i].adj[1]=0;
    487 
    488               R2 x12 = vertices[i2].r-vertices[i1].r;
    489               Real8 l12=sqrt( (x12,x12));       
    490               if (!hvertices) {
    491                 vertices[i1].color++;
    492                 vertices[i2].color++;
    493                 len[i1]+= l12;
    494                 len[i2] += l12;}
    495               hmin = Min(hmin,l12);
    496             }
    497           // definition  the default of the given mesh size
    498           if (!hvertices)
    499             {
    500             for (i=0;i<nbv;i++)
    501               if (vertices[i].color > 0)
    502                 vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
    503               else
    504                 vertices[i].m=  Metric(hmin);
    505               delete [] len;
    506             }
    507           if(verbosity>5)
    508             cout << "     hmin " << hmin << endl;
    509           // construction of edges[].adj
    510             for (i=0;i<nbv;i++)
    511               vertices[i].color = (vertices[i].color ==2) ? -1 : -2;
    512             for (i=0;i<nbe;i++)
    513               for (j=0;j<2;j++)
    514                 {
    515                   Vertex *v=edges[i].v[j];
    516                   Int4 i0=v->color,j0;
    517                   if(i0==-1)
    518                     v->color=i*2+j;
    519                   else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
    520                     j0 =  i0%2;
    521                     i0 =  i0/2;
    522                     assert( v ==  edges[i0 ].v[j0]);
    523                     edges[i ].adj[ j ] =edges +i0;
    524                     edges[i0].adj[ j0] =edges +i ;
    525                     assert(edges[i0].v[j0] == v);
    526                     //      if(verbosity>8)
    527                     //  cout << " edges adj " << i0 << " "<< j0 << " <-->  "  << i << " " << j << endl;
    528                       v->color = -3;}
    529                 }
    530 
    531         }
    532        
    533 /*   ne peut pas marche car il n'y a pas de flag dans les aretes du maillages
    534        else if (!strcmp(fieldname,"RequiredEdges"))
    535         {
    536           int i,j,n;
    537           f_in  >> n ;
    538 
    539           for (i=0;i<n;i++) {     
    540             f_in  >>  j ;
    541             assert( j <= nbe );
    542             assert( j > 0 );
    543             j--;
    544             edges[j].SetRequired();  }
    545       }
    546 */     
    547       else if (!strcmp(fieldname,"EdgeOnGeometricEdge"))
    548         {
    549           assert(edges);
    550           int i1,i2,i,j;
    551           f_in >> i2 ;
    552           if(verbosity>3)
    553             cout << "     Record EdgeOnGeometricEdge: Nb " << i2 <<endl;
    554           for (i1=0;i1<i2;i1++)
    555             {
    556               f_in  >> i >> j ;
    557               if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe))
    558                 {
    559                   cerr << "      Record EdgeOnGeometricEdge i=" << i << " j = " << j;
    560                   cerr << " nbe = " << nbe << " Gh.nbe = " <<  Gh.nbe << endl;
    561                   cerr << " We must have : (i>0 && j >0 && i <= nbe && j <= Gh.nbe) ";
    562                   cerr << " Fatal error in file " << name << " line " << f_in.LineNumber << endl;
    563                   MeshError(1);
    564                 }
    565                
    566                
    567               edges[i-1].on = Gh.edges + j-1;
    568             }
    569           //      cout << "end EdgeOnGeometricEdge" << endl;
    570         }
    571       else if (!strcmp(fieldname,"SubDomain") || !strcmp(fieldname,"SubDomainFromMesh") )
    572         {
    573          
    574           f_in >> NbSubDomains ;
    575           subdomains = new SubDomain [ NbSubDomains ];
    576             for (i=0;i< NbSubDomains;i++)
    577               { Int4 i3,head,sens;
    578               f_in  >>  i3 >>   head >> sens  >> subdomains[i].ref ;
    579                 assert (i3==3);
    580                 head --;
    581                 assert(head < nbt && head >=0);
    582                 subdomains[i].head = triangles+head;           
    583               }
    584         }
    585       else
    586         { // unkown field
    587           field = ++showfield;
    588           if(showfield==1) // just to show one time
    589             if (verbosity>5)
    590               cout << "     Warning we skip the field " << fieldname << " at line " << f_in.LineNumber << endl;
    591         }
    592       showfield=field; // just to show one time
    593     }
    594    
    595     if (ifgeom==0)
    596       {
    597         if (verbosity)
    598           cout  << " ## Warning: Mesh without geometry we construct a geometry (theta ="
    599                 << cutoffradian*180/Pi << " degres )" << endl;
    600         ConsGeometry(cutoffradian);     
    601         Gh.AfterRead();
    602       }   
    603 }
    604 
    605 
    606 void Triangles::Read_am_fmt(MeshIstream & f_in)
    607 {
    608         long int verbosity=0;
    609 
    610   Metric M1(1);
    611 
    612   if (verbosity>1)
    613     cout << "  -- ReadMesh .am_fmt file " <<  f_in.CurrentFile  << endl;
    614        
    615      Int4 i;
    616      f_in.cm() >> nbv >> nbt ;
    617      if (verbosity>3)
    618        cout << "    nbv = " << nbv  << " nbt = " << nbt << endl;
    619      f_in.eol() ;//
    620      nbvx = nbv;
    621      nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    622      triangles =new Triangle[nbtx];
    623      assert(triangles);
    624      vertices=new Vertex[nbvx];
    625      ordre=new (Vertex* [nbvx]);
    626      
    627      for (     i=0;i<nbt;i++)
    628        {
    629          Int4 i1,i2,i3;
    630          f_in >>  i1 >>  i2 >>   i3 ;
    631          triangles[i]  = Triangle(this,i1-1,i2-1,i3-1); 
    632        }
    633       f_in.eol() ;//
    634      
    635       for ( i=0;i<nbv;i++)
    636         f_in >> vertices[i].r.x >>   vertices[i].r.y,
    637           vertices[i].m = M1,vertices[i].DirOfSearch =NoDirOfSearch;
    638 
    639       f_in.eol() ;//
    640      
    641       for ( i=0;i<nbt;i++) 
    642         f_in >> triangles[i].color;
    643       f_in.eol() ;//
    644      
    645       for ( i=0;i<nbv;i++) 
    646              f_in >> vertices[i].ReferenceNumber;
    647      
    648      
    649 }
    650 
    651 ////////////////////////
    652 
    653 void  Triangles::Read_am(MeshIstream &ff)
    654 {
    655         long int verbosity=0;
    656 
    657   if (verbosity>1)
    658     cout << "  -- ReadMesh .am_fmt file " <<  ff.CurrentFile  << endl;
    659     Metric M1(1);       
    660  
    661   IFortranUnFormattedFile f_in(ff);
    662 
    663   Int4  l=f_in.Record();
    664   assert(l==2*sizeof(Int4));
    665   f_in >> nbv >> nbt ;
    666   l=f_in.Record();
    667   assert((size_t) l==nbt*sizeof(long)*4 + nbv*(2*sizeof(float)+sizeof(long)));
    668   if (verbosity>3)
    669     cout << "    nbv = " << nbv  << " nbt = " << nbt << endl;
    670  
    671   nbvx = nbv;
    672   nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    673   triangles =new Triangle[nbtx];
    674   assert(triangles);
    675   vertices=new Vertex[nbvx];
    676   ordre=new (Vertex* [nbvx]);
    677  
    678 
    679   Int4 i;
    680   for (     i=0;i<nbt;i++) {
    681     long i1,i2,i3;
    682     f_in >>  i1 >>  i2 >>   i3 ;
    683     triangles[i]  = Triangle(this,i1-1,i2-1,i3-1); }
    684  
    685   for ( i=0;i<nbv;i++) {
    686     float x,y;
    687     f_in >> x >> y;
    688     vertices[i].r.x =x;
    689     vertices[i].r.y=y;
    690     vertices[i].m=M1;}
    691  
    692   for ( i=0;i<nbt;i++) {
    693     long i;
    694     f_in >> i;
    695     triangles[i].color=i;}
    696  
    697   for ( i=0;i<nbv;i++) {
    698     long i;
    699     f_in >> i;
    700     vertices[i].ReferenceNumber=i;}
    701 }
    702 
    703 //////////////////////////////////
    704 
    705 void  Triangles::Read_nopo(MeshIstream & ff)
    706 {
    707         long int verbosity=0;
    708 
    709 
    710  if (verbosity>1)
    711     cout << "  -- ReadMesh .nopo file " <<  ff.CurrentFile  << endl;
    712  IFortranUnFormattedFile f_in(ff);
    713  
    714  
    715  Int4  l,i,j;
    716  l=f_in.Record();
    717  l=f_in.Record();
    718  f_in >> i;
    719  assert(i==32);
    720  Int4 niveau,netat,ntacm;
    721 
    722  char titre[80+1],  date[2*4+1], nomcre[6*4+1], typesd[5];
    723  f_in.read4(titre,20);
    724  f_in.read4(date,2);
    725  f_in.read4(nomcre,6);
    726  f_in.read4(typesd,1);
    727 
    728 
    729    f_in >> niveau>>netat>>ntacm;
    730  if (strcmp("NOPO",typesd))
    731    {
    732      cout << " where in record  " << f_in.where() << " " << strcmp("NOPO",typesd) << endl;
    733      cerr << " not a  nopo file but `" << typesd <<"`"<< " len = " << strlen(typesd) << endl;
    734      cerr << (int) typesd[0] << (int) typesd[1] << (int) typesd[2] << (int) typesd[3] << (int) typesd[4] << endl;
    735      cout << " nomcre :" << nomcre << endl;
    736      cout << " date   :" << date << endl;
    737      cout << " titre  :" << titre<< endl;
    738      MeshError(112);
    739    }
    740  if(verbosity>2)
    741    cout << "    nb de tableau associe : " << ntacm << " niveau =" << niveau << endl;
    742  
    743  for (i=0;i<ntacm;i++)
    744    f_in.Record();
    745 
    746  f_in.Record();
    747  f_in >> l;
    748  assert(l == 27);
    749  Int4 nop2[27];
    750  for (i=0;i<27;i++)
    751    f_in >> nop2[i];
    752  Int4 ndim = nop2[0];
    753  Int4 ncopnp = nop2[3];
    754  Int4 ne = nop2[4];
    755  Int4 ntria = nop2[7];
    756  Int4 nquad = nop2[8];
    757  Int4 np = nop2[21];
    758  // Int4 nef = nop2[13];
    759  Metric M1(1);
    760  if(verbosity>2)
    761    cout << "    ndim = " << ndim << " ncopnp= " << ncopnp << " ne = " << ne
    762         << "    ntri = " << ntria << " nquad = " << nquad << " np = " << np << endl;
    763  nbv = np;
    764  nbt = 2*nquad + ntria;
    765  if (ne != nquad+ntria || ndim !=2 || ncopnp != 1 )
    766    {
    767      cerr << " not only tria & quad in nopo mesh on dim != 2 ou ncopnp != 1 " << endl;
    768      MeshError(113);
    769    }
    770  if( nop2[24]>=0)  f_in.Record();
    771   NbOfQuad = nquad;
    772   nbvx = nbv;
    773   nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    774   triangles =new Triangle[nbtx];
    775   assert(triangles);
    776   vertices=new Vertex[nbvx];
    777   ordre=new (Vertex* [nbvx]);
    778 
    779 
    780  f_in >> l;
    781 
    782   if(verbosity>9)
    783     cout << " Read cnop4 nb of float  " << l << endl;
    784  
    785   assert(l==2*np);
    786   for (i=0;i<np;i++)
    787     {  float x,y;
    788     f_in >>  x>> y;
    789     vertices[i].r.x=x;
    790     vertices[i].r.y=y;
    791     vertices[i].m=M1;
    792     vertices[i].DirOfSearch =NoDirOfSearch;
    793 
    794     }
    795   f_in.Record();
    796   // lecture de nop5 bonjour les degats
    797   f_in >> l;
    798   if(verbosity>9)
    799     cout << " Read nop5  nb of int4 " << l << endl;
    800  Int4 k=0;
    801  Int4 nbe4 =  3*ntria + 4*nquad;
    802  // cout << " nbv = " << nbv << " nbe4 " << nbe4 << endl;
    803  SetOfEdges4 * edge4= new SetOfEdges4(nbe4,nbv);
    804  Int4 * refe = new Int4[nbe4];
    805  Int4 kr =0;
    806  for (i=0;i<ne;i++)
    807    {
    808      // Int4 ng[4]={0,0,0,0};
    809      Int4 np[4],rv[4],re[4];
    810      Int4 ncge,nmae,ndsde,npo;
    811      f_in >> ncge >> nmae >> ndsde >> npo ;
    812      //cout << " element " << i << " " << ncge << " " 
    813      // << nmae <<" " <<  npo << endl;
    814      if (ncge != 3 && ncge != 4)
    815        {
    816          cerr << " read nopo type element[" << i << "] ="
    817               << ncge << " not 3 or 4 " << endl;
    818          MeshError(115);
    819        }
    820      if (npo != 3 && npo != 4)
    821        {
    822          cerr << " read nopo element[" << i << "] npo = " 
    823               << npo << " not 3 or 4 " << endl;
    824          MeshError(115);
    825        }
    826      
    827      for( j=0;j<npo;j++)
    828        {f_in >>np[j];np[j]--;}
    829      
    830      if (ncopnp !=1)
    831        {
    832          f_in >> npo;
    833          if (npo != 3 || npo != 4)
    834            {
    835              cerr << " read nopo type element[" << i << "]= " 
    836                   << ncge << " not 3 or 4 " << endl;
    837              MeshError(115);
    838            }
    839          
    840          for(j=0;j<npo;j++)
    841            {f_in >>np[j];np[j]--;}
    842          
    843        }
    844      if (nmae>0)
    845        {
    846          Int4  ining; // no ref
    847          
    848          f_in>>ining;
    849          if (ining==1)
    850            MeshError(116);
    851          if (ining==2)
    852            for (j=0;j<npo;j++)
    853              f_in >>re[j];
    854          for (j=0;j<npo;j++)
    855            f_in >>rv[j];
    856          
    857          
    858          // set the ref on vertex and the shift np of -1 to start at 0
    859          for (j=0;j<npo;j++)
    860            vertices[np[j]].ReferenceNumber = rv[j];
    861          
    862          if (ining==2)
    863            for (j=0;j<npo;j++)
    864              if (re[j])
    865                {
    866                  kr++;
    867                  Int4 i0 = np[j];
    868                  Int4 i1 = np[(j+1)%npo];
    869                  // cout << kr << " ref  edge " << i0 << " " << i1 << " " << re[j] << endl;
    870                  refe[edge4->addtrie(i0,i1)]=re[j];
    871                }
    872        }
    873      
    874      if (npo==3)
    875        { // triangles
    876          triangles[k]  = Triangle(this,np[0],np[1],np[2]);
    877          triangles[k].color = ndsde;
    878          k++;
    879        }
    880      else if (npo==4)
    881        { // quad
    882          Triangle & t1 = triangles[k++];
    883          Triangle & t2 = triangles[k++];
    884          t1 = Triangle(this,np[0],np[1],np[2]);
    885          t2 = Triangle(this,np[2],np[3],np[0]);
    886          t1.SetHidden(OppositeEdge[1]); // two time  because the adj was not created
    887          t2.SetHidden(OppositeEdge[1]);
    888          t1.color = ndsde;
    889          t2.color = ndsde;
    890        }
    891      else
    892        {
    893          cerr << " read nopo type element =" << npo << " not 3 or 4 " << endl;
    894          MeshError(114);
    895        }
    896      
    897      
    898    }
    899  // cout << k << " == " << nbt << endl;
    900  assert(k==nbt);
    901 
    902  nbe = edge4->nb();
    903  if (nbe)
    904    {
    905      if (verbosity>7)
    906      cout << " Nb of ref edges = " << nbe << endl;
    907      if (edges)
    908        delete [] edges;
    909      edges = new Edge[nbe];
    910      for (i=0;i<nbe;i++)
    911        {
    912          edges[i].v[0] = vertices + edge4->i(i);
    913          edges[i].v[1] = vertices + edge4->j(i);
    914          edges[i].ref = refe[i];
    915          //      cout << i << " " <<  edge4->i(i) << " " <<  edge4->j(i) << endl;
    916        }
    917       if (verbosity>7)
    918         cout << " Number of reference edge in the  mesh = " << nbe << endl;
    919    }
    920  delete [] refe;
    921  delete edge4;
    922 }
    923   void  Triangles::Read_ftq(MeshIstream & f_in)
    924 {
    925         long int verbosity=0;
    926 
    927   // 
    928   if (verbosity>1)
    929     cout << "  -- ReadMesh .ftq file " <<  f_in.CurrentFile  << endl;
    930  
    931   Int4 i,ne,nt,nq;
    932   f_in.cm() >> nbv >> ne >> nt >> nq ;
    933   if (verbosity>3)
    934     cout << "    nbv = " << nbv  << " nbtra = " << nt << " nbquad = " << nq << endl;
    935   nbt = nt+2*nq;
    936  
    937 
    938   nbvx = nbv;
    939   nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    940   triangles =new Triangle[nbtx];
    941   assert(triangles);
    942   vertices=new Vertex[nbvx];
    943   ordre=new (Vertex* [nbvx]);
    944   Int4 k=0;
    945  
    946   for ( i=0;i<ne;i++)
    947     {
    948       long ii,i1,i2,i3,i4,ref;
    949       f_in >>  ii;
    950         if (ii==3)
    951           { // triangles
    952             f_in >> i1>>  i2 >>   i3 >> ref ;
    953             triangles[k]  = Triangle(this,i1-1,i2-1,i3-1);
    954             triangles[k++].color = ref;
    955           }
    956         else if (ii==4)
    957           { // quad
    958             f_in >> i1>>  i2 >>   i3 >> i4 >> ref ;
    959             Triangle & t1 = triangles[k++];
    960             Triangle & t2 = triangles[k++];
    961             t1 = Triangle(this,i1-1,i2-1,i3-1);
    962             t1.color=ref;
    963             t2 = Triangle(this,i3-1,i4-1,i1-1);
    964             t2.color=ref;
    965             t1.SetHidden(OppositeEdge[1]); // two time  because the adj was not created
    966             t2.SetHidden(OppositeEdge[1]);       
    967           }
    968         else
    969           {
    970             cout << " read ftq type element =" << ii << " not 3 or 4 " << endl;
    971             MeshError(111);
    972           }
    973     }
    974   assert(k==nbt);
    975   Metric M1(1);
    976   for ( i=0;i<nbv;i++)
    977     {
    978       f_in >> vertices[i].r.x >>   vertices[i].r.y >> vertices[i].ReferenceNumber;
    979        vertices[i].DirOfSearch =NoDirOfSearch;
    980        vertices[i].m = M1;
    981 
    982     }
    983 }
    984 
    985 ///////////////////////////////////////////////
    986 
    987 void  Triangles::Read_msh(MeshIstream &f_in)
    988 {
    989         long int verbosity=0;
    990 
    991     Metric M1(1.);
    992   if (verbosity>1)
    993     cout << "  -- ReadMesh .msh file " <<  f_in.CurrentFile  << endl;
    994        
    995      Int4 i;
    996      f_in.cm() >> nbv >> nbt ;
    997      while (f_in.in.peek()==' ')
    998          f_in.in.get();
    999      if(isdigit(f_in.in.peek()))
    1000        f_in >> nbe;
    1001      if (verbosity>3)
    1002        cout << "    nbv = " << nbv  << " nbt = " << nbt << " nbe = " << nbe << endl;
    1003      nbvx = nbv;
    1004      nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    1005      triangles =new Triangle[nbtx];
    1006      assert(triangles);
    1007      vertices=new Vertex[nbvx];
    1008      ordre=new (Vertex* [nbvx]);
    1009       edges = new Edge[nbe];
    1010      for ( i=0;i<nbv;i++)
    1011         {
    1012          f_in >> vertices[i].r.x >>   vertices[i].r.y >> vertices[i].ReferenceNumber;
    1013             vertices[i].on=0;
    1014             vertices[i].m=M1;
    1015          //if(vertices[i].ReferenceNumber>NbRef)        NbRef=vertices[i].ReferenceNumber; 
    1016         }
    1017      for (     i=0;i<nbt;i++)
    1018        {
    1019          Int4 i1,i2,i3,r;
    1020          f_in >>  i1 >>  i2 >>   i3 >> r;
    1021          triangles[i]  = Triangle(this,i1-1,i2-1,i3-1);
    1022          triangles[i].color = r;         
    1023        }
    1024      for (i=0;i<nbe;i++)
    1025        {
    1026          Int4 i1,i2,r;
    1027          f_in >>  i1 >>  i2 >> r;
    1028               edges[i].v[0]= vertices +i1-1;
    1029               edges[i].v[1]= vertices +i2-1;
    1030               edges[i].adj[0]=0;
    1031               edges[i].adj[1]=0;
    1032               edges[i].ref = r;
    1033               edges[i].on=0;
    1034        }
    1035      
    1036 }
    1037 
    1038 //////////////////////////////////////////////////
    1039 
    1040 void  Triangles::Read_amdba(MeshIstream &f_in )
    1041 {
    1042         long int verbosity=0;
    1043 
    1044   Metric M1(1);
    1045   if (verbosity>1)
    1046     cout << "  -- ReadMesh .amdba file " <<  f_in.CurrentFile  << endl;
    1047        
    1048      Int4 i;
    1049      f_in.cm() >> nbv >> nbt ;
    1050      //    if (verbosity>3)
    1051        cout << "    nbv = " << nbv  << " nbt = " << nbt << endl;
    1052      f_in.eol() ;//
    1053      nbvx = nbv;
    1054      nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals
    1055      triangles =new Triangle[nbtx];
    1056      assert(triangles);
    1057      vertices=new Vertex[nbvx];
    1058      ordre=new (Vertex* [nbvx]);
    1059      Int4 j;
    1060       for ( i=0;i<nbv;i++)
    1061         {
    1062           f_in >> j ;
    1063           assert( j >0 && j <= nbv);
    1064           j--;
    1065           f_in >> vertices[j].r.x >>   vertices[j].r.y >> vertices[j].ReferenceNumber;
    1066            vertices[j].m=M1;
    1067            vertices[j].DirOfSearch=NoDirOfSearch;
    1068         }
    1069      
    1070       for (     i=0;i<nbt;i++)
    1071        {
    1072          Int4 i1,i2,i3,ref;
    1073            f_in >> j ;
    1074            assert( j >0 && j <= nbt);
    1075            j--;
    1076            f_in >> i1 >>  i2 >>   i3 >> ref;
    1077          triangles[j]  = Triangle(this,i1-1,i2-1,i3-1);
    1078          triangles[j].color =ref;
    1079        }
    1080       f_in.eol() ;//
    1081      
    1082   // cerr << " a faire " << endl;
    1083   //MeshError(888);
    1084 }
    1085 
    1086 
    1087 Triangles::Triangles(const char * filename,Real8 cutoffradian)
    1088 : Gh(*(new Geometry())), BTh(*this)
    1089 {
    1090 
    1091   //  Int4 beginquad=0,begintria=0;
    1092   // Int4 endquad=0;endtria=0;
    1093   //int type_file=0;
    1094 
    1095   int lll = strlen(filename);
    1096   int  am_fmt = !strcmp(filename + lll - 7,".am_fmt");
    1097   int  amdba = !strcmp(filename + lll - 6,".amdba");
    1098   int  am = !strcmp(filename + lll - 3,".am");
    1099   int  nopo = !strcmp(filename + lll - 5,".nopo");
    1100   int  msh = !strcmp(filename + lll - 4,".msh");
    1101   int  ftq = !strcmp(filename + lll - 4,".ftq");
    1102 
    1103  // cout << " Lecture type  :" << filename + lll - 7 <<":" <<am_fmt<<  endl;
    1104 
    1105   char * cname = new char[lll+1];
    1106   strcpy(cname,filename);
    1107   Int4 inbvx =0;
    1108   PreInit(inbvx,cname);
    1109   OnDisk = 1;
    1110   //  allocGeometry = &Gh; // after Preinit ;
    1111 
    1112   MeshIstream f_in (filename);
    1113  
    1114   if (f_in.IsString("MeshVersionFormatted"))
    1115     {
    1116       int version ;
    1117       f_in >> version ;
    1118       Read(f_in,version,cutoffradian);
    1119     }
    1120   else {     
    1121     if (am_fmt) Read_am_fmt(f_in);
    1122     else if (am) Read_am(f_in);
    1123     else if (amdba) Read_amdba(f_in);
    1124     else if (msh) Read_msh(f_in);
    1125     else if (nopo) Read_nopo(f_in);
    1126     else if (ftq) Read_ftq(f_in);
    1127     else
    1128       {
    1129         cerr << " Unkown type mesh " << filename << endl;
    1130         MeshError(2);
    1131       }
    1132       ConsGeometry(cutoffradian);
    1133       Gh.AfterRead();   
    1134   }
    1135 
    1136   SetIntCoor();
    1137   FillHoleInMesh();
    1138   // Make the quad ---
    1139    
    1140 }
    1141 
    1142270Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){
    1143271
    … …  
    1148276  SetIntCoor();
    1149277  FillHoleInMesh();
    1150 }
    1151 
    1152 void Geometry::ReadGeometry(const char * filename) {
    1153         long int verbosity=0;
    1154 
    1155   OnDisk = 1;
    1156   if(verbosity>1)
    1157       cout << "  -- ReadGeometry " << filename << endl;
    1158   MeshIstream f_in (filename);
    1159   ReadGeometry(f_in,filename);
    1160278}
    1161279

  • issm/trunk/src/c/Bamgx/MeshWrite.cpp

    r2794 r2795  
    1414
    1515namespace bamg {
    16 
    17 void Triangles::Write(const char * filename,const TypeFileMesh typein )
    18 {
    19         long int verbosity=0;
    20 
    21   TypeFileMesh type = typein;
    22   const char * gsuffix=".gmsh";
    23   int ls=0;
    24   int lll = strlen(filename);
    25   if (type==AutoMesh)
    26     {
    27       type = BDMesh;
    28       if      (!strcmp(filename + lll - (ls=7),".am_fmt")) type = am_fmtMesh;
    29       else if (!strcmp(filename + lll - (ls=6),".amdba"))  type = amdbaMesh;
    30       else if (!strcmp(filename + lll - (ls=3),".am"))     type = amMesh;
    31       else if (!strcmp(filename + lll - (ls=5),".nopo"))   type = NOPOMesh;
    32       else if (!strcmp(filename + lll - (ls=4),".msh"))    type = mshMesh;
    33       else if (!strcmp(filename + lll - (ls=4),".ftq"))    type = ftqMesh;
    34       else if (!strcmp(filename + lll - (ls=7),".AM_FMT")) type = am_fmtMesh;
    35       else if (!strcmp(filename + lll - (ls=6),".AMDBA"))  type = amdbaMesh;
    36       else if (!strcmp(filename + lll - (ls=3),".AM"))     type = amMesh;
    37       else if (!strcmp(filename + lll - (ls=5),".NOPO"))   type = NOPOMesh;
    38       else if (!strcmp(filename + lll - (ls=4),".MSH"))    type = mshMesh;
    39       else if (!strcmp(filename + lll - (ls=4),".FTQ"))    type = ftqMesh;
    40       else ls=0;
    41     }
    42   if (verbosity>1)
    43     {
    44       cout << "  -- Writing the file " << filename << " of type " ;
    45      switch (type)
    46        {
    47        case BDMesh     :  cout << " BD Mesh "  ; break;
    48        case NOPOMesh   :  cout << " NOPO "     ; break;
    49        case amMesh     :  cout << " am "       ; break;
    50        case am_fmtMesh :  cout << " am_fmt "   ; break;
    51        case amdbaMesh  :  cout << " amdba "    ; break;
    52        case ftqMesh    :  cout << " ftq "      ; break;
    53        case mshMesh    :  cout << " msh "      ; break;
    54         default:
    55           cerr << endl
    56                <<  " Unkown type mesh file " << (int) type << " for Writing " << filename <<endl;
    57           MeshError(1);
    58        }     
    59      Int4 NbOfTria =  nbt-2*NbOfQuad-NbOutT ;
    60      if (NbOfTria)      cout << " NbOfTria = " << NbOfTria;
    61      if (NbOfQuad)      cout << " NbOfQuad = " << NbOfQuad;
    62      if (nbe)           cout << " NbOfRefEdge = " << nbe ;
    63      cout    << endl;
    64 
    65     }
    66   ofstream f(filename /*,ios::trunc*/);
    67  f.precision(12);
    68 
    69    if (f)
    70      switch (type)
    71        {
    72        case BDMesh     :
    73           {
    74              if ( ! Gh.OnDisk)
    75               {
    76                  delete [] Gh.name;
    77                  Gh.name = new char[lll+1+strlen(gsuffix)];
    78                  strcpy(Gh.name,filename);
    79                  if (Gh.name[lll-ls-1]=='.') strcpy(Gh.name+lll-ls, gsuffix+1);
    80                  else strcpy(Gh.name+lll-ls,gsuffix);
    81                 cout << " write geo in " << Gh.name << endl;
    82                   ofstream f(Gh.name) ;
    83                   f << Gh ;
    84                   Gh.OnDisk=true;
    85               }
    86              f << *this     ;
    87              break;
    88            }
    89        case NOPOMesh   :  Write_nopo(f)  ; break;
    90        case amMesh     :  Write_am(f)    ; break;
    91        case am_fmtMesh :  Write_am_fmt(f); break;
    92        case amdbaMesh  :  Write_amdba(f) ; break;
    93        case ftqMesh    :  Write_ftq(f)   ; break;
    94        case mshMesh    :  Write_msh(f)   ; break;
    95         default:
    96           cerr << " Unkown type mesh file " << (int) type << " for Writing " << filename <<endl;
    97           MeshError(1);
    98        }
    99    else
    100      {
    101        cerr << " Error openning file " << filename << endl;
    102        MeshError(1);
    103      }
    104    if(verbosity>5)
    105    cout << "end write" << endl;
    106      
    107 }
    108 void Triangles::Write_nop5(OFortranUnFormattedFile * f,
    109                            Int4 &lnop5,Int4 &nef,Int4 &lgpdn,Int4 ndsr) const
    110 {
    111         long int verbosity=0;
    112 
    113   ndsr =0;
    114   Int4 * reft = new Int4[nbt];
    115   //Int4 nbInT = ;
    116   ConsRefTriangle(reft);
    117   Int4 no5l[20];
    118 
    119   Int4 i5 = 0;
    120   Int4 i,j,k=0,l5;
    121   //  Int4 ining=0;
    122   Int4 imax,imin;
    123 
    124   lgpdn = 0;
    125   nef=0;
    126   // construction of a liste linked  of edge
    127   Edge ** head  = new Edge *[nbv];
    128   Edge  ** link = new Edge * [nbe];
    129   for (i=0;i<nbv;i++)
    130     head[i]=0; // empty liste
    131  
    132   for (i=0;i<nbe;i++)
    133     {
    134       j = Min(Number(edges[i][0]),Number(edges[i][1]));
    135       link[i]=head[j];
    136       head[j]=edges +i;
    137     }
    138   for ( i=0;i<nbt;i++)
    139     {
    140       no5l[0]    = 0;
    141       Int4 kining=0;
    142       Int4 ining=0;
    143       Int4 nmae =0;
    144       Int4 np=0;
    145       l5 = 2;
    146       Triangle & t = triangles[i];
    147       Triangle * ta; //
    148       Vertex *v0,*v1,*v2,*v3;
    149       if (reft[i]<0) continue;
    150       ta = t.Quadrangle(v0,v1,v2,v3);
    151       if (!ta)
    152         { // a triangles
    153           no5l[l5++] = Max(subdomains[reft[i]].ref,(Int4) 1);
    154           np = 3;
    155           no5l[l5++] = np;
    156           no5l[0]    = np;
    157           no5l[l5++] = Number(triangles[i][0]) +1;
    158           no5l[l5++] = Number(triangles[i][1]) +1;
    159           no5l[l5++] = Number(triangles[i][2]) +1;
    160           imax = Max3(no5l[l5-1],no5l[l5-2],no5l[l5-3]);
    161           imin = Min3(no5l[l5-1],no5l[l5-2],no5l[l5-3]);
    162           lgpdn = Max(lgpdn,imax-imin);
    163           kining=l5++;
    164           // ref of 3 edges
    165           for (j=0;j<3;j++)
    166             {
    167               no5l[l5] = 0;
    168               int i0 = (int) j;
    169               int i1 = (i0+1) %3;
    170               Int4 j1= no5l[4+i0];
    171               Int4 j2= no5l[4+i1];
    172               Int4 ji = Min(j1,j2)-1;
    173               Int4 ja = j1+j2-2;
    174               Edge * e=head[ji];
    175               while (e)
    176                 if(Number((*e)[0])+Number((*e)[1]) == ja)
    177                    {
    178                      no5l[l5] = e->ref;
    179                      break;
    180                    }
    181                 else
    182                   e = link[Number(e)];
    183               l5++;
    184             }
    185           if ( no5l[l5-1] || no5l[l5-2] || no5l[l5-3]  )
    186             ining=2;
    187           else
    188             l5 -= 3;
    189          
    190           no5l[l5++] = triangles[i][0].ref();
    191           no5l[l5++] = triangles[i][1].ref();
    192           no5l[l5++] = triangles[i][2].ref();
    193           if (ining ||  no5l[l5-1] || no5l[l5-2] || no5l[l5-3]  )
    194             ining= ining ? ining : 3;
    195           else
    196             l5 -= 3,ining=0;
    197 
    198          
    199         }
    200       else if ( &t<ta)
    201         {
    202           k++;
    203           no5l[l5++] = Max(subdomains[reft[i]].ref,(Int4) 1);
    204           np =4;
    205           no5l[l5++] = np;
    206           no5l[0]    = np;
    207          
    208           no5l[l5++] = Number(v0) +1;
    209           no5l[l5++] = Number(v1) +1;
    210           no5l[l5++] = Number(v2) +1;
    211           no5l[l5++] = Number(v3) +1;
    212          
    213           imax = Max(Max(no5l[l5-1],no5l[l5-2]),Max(no5l[l5-3],no5l[l5-4]));
    214           imin = Min(Min(no5l[l5-1],no5l[l5-2]),Min(no5l[l5-3],no5l[l5-4]));
    215           lgpdn = Max(lgpdn,imax-imin);
    216 
    217 
    218           kining=l5++;
    219           // ref of the 4 edges
    220           // ref of 3 edges
    221           for (j=0;j<4;j++)
    222             {
    223               no5l[l5] = 0;
    224               int i0 = (int) j;
    225               int i1 = (i0+1) %4;
    226               Int4 j1= no5l[4+i0];
    227               Int4 j2= no5l[4+i1];
    228               Int4 ji = Min(j1,j2)-1;
    229               Int4 ja = j1+j2-2;
    230               Edge *e=head[ji];
    231               while (e)
    232                 if(Number((*e)[0])+Number((*e)[1]) == ja)
    233                    {
    234                      no5l[l5] = e->ref;
    235                      break;
    236                    }
    237                 else
    238                   e = link[Number(e)];
    239               l5++;
    240             }
    241           if ( no5l[l5-1] || no5l[l5-2] || no5l[l5-3] || no5l[l5-4] )
    242             ining=2;
    243           else
    244             l5 -= 4;
    245 
    246           no5l[l5++] = v0->ref();
    247           no5l[l5++] = v1->ref();
    248           no5l[l5++] = v2->ref();
    249           no5l[l5++] = v3->ref();
    250           if (ining || no5l[l5-1] || no5l[l5-2] || no5l[l5-3] || no5l[l5-4] )
    251             ining= ining ? ining : 3;
    252           else
    253             l5 -= 4;
    254 
    255         }
    256       else l5=0;
    257 
    258      if (l5)
    259        {
    260          if (ining)
    261            {
    262              nef ++;
    263              nmae = l5-kining;
    264              no5l[kining] = ining;
    265            }
    266          else l5--;
    267          no5l[1]=nmae;
    268          // for all ref 
    269          for (j=kining+1;j<l5;j++)
    270            {
    271              no5l[j] = Abs(no5l[j]);
    272              ndsr = Max(ndsr,no5l[j]);
    273            }
    274          
    275          if (f && i < 10 && verbosity > 10)
    276            {
    277              cout << " e[ " << i << " @" << i5 << "]=";
    278              for (j=0;j<l5;j++)
    279                cout << " " << no5l[j];
    280              cout << endl;
    281            }
    282          
    283          if (f)
    284            for (j=0;j<l5;j++)
    285              *f << no5l[j];
    286          i5 += l5;
    287        }
    288     }
    289   if(verbosity>10)
    290   cout << "   fin write nopo 5 i5=" << i5 << " " << i5*4 << endl;
    291   lnop5=i5;
    292   lgpdn++; // add 1
    293   delete [] reft;
    294   delete [] head;
    295   delete [] link;
    296  
    297 }
    298 
    299 void Triangles::Write_nopo(ostream &ff) const
    300 
    301 {
    302  Int4  nef=0;
    303  Int4 lpgdn=0;
    304  Int4 ndsr=0;
    305  Int4 i;
    306  Int4 ndsd=1;
    307  Int4 lnop5=0;
    308  
    309  OFortranUnFormattedFile f(ff);
    310  
    311  for (i=0;i<NbSubDomains ;i++)
    312    ndsd=Max(ndsd,subdomains[i].ref);
    313  
    314  // to compute the lnop5,nef,lpgdn,ndsr parameter
    315  Write_nop5(0,lnop5,nef,lpgdn,ndsr);
    316  
    317  f.Record();
    318  
    319  f <<  Int4(13)<<Int4(6)<<Int4(32)<<Int4(0)<<Int4(27)<<Int4(0) ;
    320  f << Int4(nbv+nbv) ;
    321  f << lnop5;
    322  f << Int4(1 )<<Int4(1)<<Int4(1 )<<Int4(1)<<Int4(2)<<Int4(1);
    323 
    324  f.Record(33*sizeof(Int4));
    325 
    326  f << Int4(32) ;
    327  
    328  //char *c=identity;
    329  time_t timer =time(0);
    330  char buf[10];
    331  strftime(buf ,10,"%y/%m/%d",localtime(&timer));
    332  f.write4(identity,20);
    333  f.write4(buf,2);
    334  f.write4("created with BAMG",6);
    335  f.write4("NOPO",1);
    336 
    337  
    338  f << Int4(0) << Int4(1) << Int4(0) ;
    339  f.Record();
    340  Int4 nbquad= NbOfQuad;
    341  Int4 nbtria= nbt-NbOutT - 2*NbOfQuad;
    342 
    343  cout << " lnop5      = " << lnop5 << endl;
    344  cout << " nbquad     = " << nbquad << endl;
    345  cout << " nbtrai     = " << nbtria << endl;
    346  cout << " lpgdn      = " << lpgdn << endl;
    347  cout << " nef        = " << nef  << endl;
    348  cout << " np         = " << nbv  << endl;
    349  cout << " ndsr       = " << ndsr << endl;
    350  f << Int4(27) 
    351    << Int4(2)  << ndsr     << ndsd    << Int4(1) << nbtria+nbquad
    352    << Int4(0)  << Int4(0)  << nbtria  << nbquad  << Int4(0)
    353    << Int4(0)  << Int4(0)  << Int4(0) << nef     << Int4(nbv)
    354    << Int4(0)  << Int4(0)  << Int4(0) << Int4(0) << Int4(0)
    355    << Int4(0)  << nbv      << Int4(2) << lpgdn   << Int4(0)
    356    << lnop5    << Int4(1) ;
    357   f.Record();
    358   f << (Int4) 2*nbv;
    359   for (i=0;i<nbv;i++)
    360     f << (float) vertices[i].r.x <<  (float) vertices[i].r.y;
    361   f.Record();
    362   f << lnop5;
    363   Write_nop5(&f,lnop5,nef,lpgdn,ndsr);
    364   // cout << "fin write nopo" << endl;
    365 }
    366 
    367 void Triangles::Write_am_fmt(ostream &f) const
    368 {
    369   Int4 i,j;
    370   assert(this && nbt);
    371   Int4 * reft = new Int4[nbt];
    372   Int4 nbInT =    ConsRefTriangle(reft);
    373   f.precision(12);
    374   f << nbv << " " << nbInT << endl;
    375   for (i=0;i<nbt;i++)
    376     if(reft[i]>=0)
    377       {
    378         f << Number(triangles[i][0]) +1 << " " ;
    379         f << Number(triangles[i][1]) +1 << " " ;
    380         f << Number(triangles[i][2]) +1 << " " ;
    381         f << endl;
    382       }
    383   for (i=0;i<nbv;i++)
    384       f << vertices[i].r.x << " " << vertices[i].r.y << endl;
    385    for (j=i=0;i<nbt;i++)
    386      if (reft[i]>=0)
    387        f << subdomains[reft[i]].ref  << (j++%10 == 9 ?  '\n' : ' ');
    388    f << endl;
    389    for (i=0;i<nbv;i++)
    390      f << vertices[i].ref()  << (i%10 == 9 ?  '\n' : ' ');
    391    f << endl;
    392    delete [] reft;
    393 
    394 
    395 }
    396 
    397 void Triangles::Write_am(ostream &ff) const
    398 {
    399   OFortranUnFormattedFile f(ff); 
    400   Int4 i,j;
    401   assert(this && nbt);
    402   Int4 * reft = new Int4[nbt];
    403   Int4 nbInT =    ConsRefTriangle(reft);
    404   f.Record();
    405   f << nbv << nbInT ;
    406   f.Record();
    407   for (i=0;i<nbt;i++)
    408     if(reft[i]>=0)
    409       {
    410         f << Number(triangles[i][0]) +1 ;
    411         f << Number(triangles[i][1]) +1 ;
    412         f << Number(triangles[i][2]) +1 ;
    413       }
    414   for (i=0;i<nbv;i++)
    415     {
    416       float x= vertices[i].r.x;
    417       float y= vertices[i].r.y;
    418       f << x << y ;
    419     }
    420   for (j=i=0;i<nbt;i++)
    421     if (reft[i]>=0)
    422       f << subdomains[reft[i]].ref;
    423   for (i=0;i<nbv;i++)
    424     f << vertices[i].ref() ;
    425   delete [] reft;
    426 }
    427 
    428 void Triangles::Write_ftq(ostream &f) const
    429 {
    430 
    431   Int4 i;
    432   assert(this && nbt);
    433   Int4 * reft = new Int4[nbt];
    434   Int4 nbInT =    ConsRefTriangle(reft);
    435   f.precision(12);
    436   Int4 nele = nbInT-NbOfQuad;
    437   Int4 ntri =  nbInT-2*NbOfQuad;
    438   Int4 nqua =  NbOfQuad;
    439 
    440   f << nbv << " " << nele << " " << ntri <<  " " << nqua << endl;
    441   Int4 k=0;
    442   for( i=0;i<nbt;i++)
    443     {
    444       Triangle & t = triangles[i];
    445       Triangle * ta; //
    446       Vertex *v0,*v1,*v2,*v3;
    447       if (reft[i]<0) continue;
    448       ta = t.Quadrangle(v0,v1,v2,v3);
    449       if (!ta)
    450         { // a triangles
    451           f << "3 "
    452             << Number(triangles[i][0]) +1 << " "
    453             << Number(triangles[i][1]) +1 << " "
    454             << Number(triangles[i][2]) +1 << " "
    455             << subdomains[reft[i]].ref << endl;
    456           k++;
    457         }
    458       if ( &t<ta)
    459         {
    460           k++;
    461           f << "4 " << Number(v0)+1 << " " << Number(v1)+1  << " " 
    462             << Number(v2)+1 << " "  << Number(v3)+1 << " " 
    463             << subdomains[reft[i]].ref << endl;
    464         }
    465     }
    466   assert(k == nele);
    467  
    468   for (i=0;i<nbv;i++)
    469     f << vertices[i].r.x << " " << vertices[i].r.y
    470       << " " <<  vertices[i].ref() << endl;
    471   delete [] reft;
    472  
    473  
    474 }
    475 void Triangles::Write_msh(ostream &f) const
    476 {
    477   Int4 i;
    478   assert(this && nbt);
    479   Int4 * reft = new Int4[nbt];
    480   Int4 nbInT =    ConsRefTriangle(reft);
    481   f.precision(12);
    482   f << nbv << " " << nbInT << " " << nbe <<  endl;
    483 
    484   for (i=0;i<nbv;i++)
    485     f << vertices[i].r.x << " " << vertices[i].r.y << " "
    486       << vertices[i].ref() <<   endl;
    487 
    488   for (i=0;i<nbt;i++)
    489     if(reft[i]>=0)
    490       f << Number(triangles[i][0]) +1 << " "
    491         << Number(triangles[i][1]) +1 << " "
    492         << Number(triangles[i][2]) +1 << " "
    493         << subdomains[reft[i]].ref << endl;
    494  
    495 
    496   for (i=0;i<nbe;i++)
    497     f << Number(edges[i][0]) +1 << " "  << Number(edges[i][1]) +1
    498       << " " << edges[i].ref << endl;
    499      
    500    delete [] reft;
    501 
    502 }
    503 
    504 void Triangles::Write_amdba(ostream &f) const
    505 {
    506   assert(this && nbt);
    507 
    508   Int4 i,j;
    509   Int4 * reft = new Int4[nbt];
    510   Int4 nbInT =    ConsRefTriangle(reft);
    511   f << nbv << " " << nbInT << endl;
    512   cout.precision(12);
    513   for (i=0;i<nbv;i++)
    514     f << i+1 << " "
    515       << vertices[i].r.x
    516       << " " << vertices[i].r.y
    517       << " " << vertices[i].ref() << endl;
    518   j=1;
    519   for (i=0;i<nbt;i++)
    520     if(reft[i]>=0)
    521         f << j++ << " "
    522           << Number(triangles[i][0]) +1 << " "
    523           << Number(triangles[i][1]) +1 << " "
    524           << Number(triangles[i][2]) +1 << " "
    525           << subdomains[reft[i]].ref  << endl ;
    526         f << endl;
    527    delete [] reft;
    528 
    529 
    530 }
    53116
    53217void Triangles::WriteBamgMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
    … …  
    729214        }
    730215}
    731 
    732 void Triangles::Write(const char * filename)
    733 {
    734   ofstream f(filename);
    735   if (f)
    736     {
    737        if (name) delete name;
    738        name = new char[strlen(filename)+1];
    739        strcpy(name,filename);
    740        OnDisk =1;
    741        f << *this;
    742     }
    743 }
    744 void Triangles::WriteElements(ostream& f,Int4 * reft ,Int4 nbInT) const
    745    {
    746      const Triangles & Th= *this;
    747          long int verbosity=0;
    748 
    749      // do triangle and quad
    750      if(verbosity>9)
    751        cout  << " In Triangles::WriteElements " << endl
    752              << "   Nb of In triangles " << nbInT-Th.NbOfQuad*2 << endl
    753              << "   Nb of Quadrilaterals " <<  Th.NbOfQuad << endl
    754              << "   Nb of in+out+quad  triangles " << Th.nbt << " " << nbInT << endl;
    755          
    756      Int4 k=nbInT-Th.NbOfQuad*2;
    757      Int4 num =0;
    758      if (k>0) {
    759        f << "\nTriangles\n"<< k << endl;
    760        for(Int4 i=0;i<Th.nbt;i++)
    761          {
    762            Triangle & t = Th.triangles[i];
    763            if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) ))
    764             { k--;
    765               f << Th.Number(t[0])+1 << " " << Th.Number(t[1])+1
    766                    << " "  << Th.Number(t[2])+1  << " " << Th.subdomains[reft[i]].ref << endl;
    767                    reft[i] = ++num;
    768                  }
    769          }
    770      }
    771      if (Th.NbOfQuad>0) {
    772        f << "\nQuadrilaterals\n"<<Th.NbOfQuad << endl;
    773        k = Th.NbOfQuad;
    774        for(Int4 i=0;i<Th.nbt;i++)
    775          {
    776            Triangle & t = Th.triangles[i];
    777            Triangle * ta; //
    778            Vertex *v0,*v1,*v2,*v3;
    779            if (reft[i]<0) continue;
    780            if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta)
    781               {
    782                 k--;
    783                 f << Th.Number(v0)+1 << " " << Th.Number(v1)+1  << " " 
    784                   << Th.Number(v2)+1 << " "  << Th.Number(v3)+1 << " " 
    785                   << Th.subdomains[reft[i]].ref << endl;
    786                   reft[i] = ++num;
    787                   reft[Number(ta)] = num;
    788               }
    789          }
    790        assert(k==0);
    791      }
    792      // warning reft is now the element number
    793    }
    794 
    795 ostream& operator <<(ostream& f, const   Triangles & Th)
    796  {
    797   //  Th.FindSubDomain();
    798    // warning just on say the class is on the disk
    799   //  ((Triangles *) &Th)->OnDisk = 1;
    800 
    801    Int4 * reft = new Int4[Th.nbt];
    802    Int4 nbInT =    Th.ConsRefTriangle(reft);
    803    {
    804      f << "MeshVersionFormatted 0" <<endl;
    805      f << "\nDimension\n"  << 2 << endl;
    806      f << "\nIdentifier\n" ;
    807      WriteStr(f,Th.identity);
    808      f << "\n\nGeometry\n" ;
    809      if( Th.Gh.OnDisk)
    810        WriteStr(f,Th.Gh.name),     f <<endl;
    811      else
    812        { // empty file name -> geom in same file
    813          f << "\"\"" << endl << endl;
    814          f << "# BEGIN of the include geometry file because geometry is not on the disk"
    815            << Th.Gh << endl;
    816          f << "End" << endl
    817            << "# END of the include geometrie file because geometry is not on the disk"
    818            << endl ;
    819        }
    820    }
    821    {
    822      f.precision(12);
    823      f << "\nVertices\n" << Th.nbv <<endl;
    824      for (Int4 i=0;i<Th.nbv;i++)
    825        {
    826          Vertex & v =  Th.vertices[i];
    827          f << v.r.x << " " << v.r.y << " " << v.ref() << endl;
    828        }
    829    }
    830   Int4 ie;
    831    {
    832      f << "\nEdges\n"<< Th.nbe << endl;
    833      for(ie=0;ie<Th.nbe;ie++)
    834        {
    835          Edge & e = Th.edges[ie];
    836          f << Th.Number(e[0])+1 << " " << Th.Number(e[1])+1;
    837         f << " " << e.ref <<endl;
    838        }
    839      if(Th.NbCrackedEdges)
    840        {
    841          f << "\nCrackedEdges\n"<< Th.NbCrackedEdges << endl;
    842          for( ie=0;ie<Th.NbCrackedEdges;ie++)
    843            {
    844              Edge & e1 = *Th.CrackedEdges[ie].a.edge;
    845              Edge & e2 = *Th.CrackedEdges[ie].b.edge;
    846              f << Th.Number(e1)+1 << " " << Th.Number(e2)+1 <<endl;;
    847            }
    848        }
    849    }
    850 
    851    Th.WriteElements(f,reft,nbInT);
    852    {
    853      f << "\nSubDomainFromMesh\n" << Th.NbSubDomains<< endl ;
    854      for (Int4 i=0;i<Th.NbSubDomains;i++)
    855        f << 3 << " " << reft[Th.Number(Th.subdomains[i].head)] << " " << 1 << " "
    856          <<  Th.subdomains[i].ref << endl;
    857      
    858    }
    859    if (Th.Gh.NbSubDomains)
    860      {
    861         f << "\nSubDomainFromGeom\n" << Th.Gh.NbSubDomains << endl ;
    862        for (Int4 i=0;i<Th.NbSubDomains;i++)
    863          { 
    864          f << 2 << " " << Th.Number(Th.subdomains[i].edge)+1 << " "
    865            <<  Th.subdomains[i].sens  << " " <<  Th.Gh.subdomains[i].ref << endl;
    866          }
    867    }
    868    {
    869      f << "\nVertexOnGeometricVertex\n"<<  Th.NbVerticesOnGeomVertex << endl;
    870      for (Int4 i0=0;i0<Th.NbVerticesOnGeomVertex;i0++)
    871        {
    872          VertexOnGeom & v =Th.VerticesOnGeomVertex[i0];
    873          assert(v.OnGeomVertex()) ;
    874          f << " " << Th.Number(( Vertex *)v)+1 
    875            << " " << Th.Gh.Number(( GeometricalVertex * )v)+1
    876            << endl;
    877        }
    878    }
    879    {
    880      f << "\nVertexOnGeometricEdge\n"<<  Th.NbVerticesOnGeomEdge << endl;
    881      for (Int4 i0=0;i0<Th.NbVerticesOnGeomEdge;i0++)
    882        {
    883          const VertexOnGeom & v =Th.VerticesOnGeomEdge[i0];
    884          assert(v.OnGeomEdge()) ;   
    885          f << " " << Th.Number((Vertex * )v)+1  ;
    886          f << " " << Th.Gh.Number((const  GeometricalEdge * )v)+1  ;
    887          f << " " << (Real8 ) v << endl;
    888        }
    889    }
    890    {
    891      Int4 i0,k=0;
    892 
    893      for (i0=0;i0<Th.nbe;i0++)
    894        if ( Th.edges[i0].on ) k++;
    895      
    896      f << "\nEdgeOnGeometricEdge\n"<< k << endl;
    897       for (i0=0;i0<Th.nbe;i0++)
    898                  if ( Th.edges[i0].on )
    899                   f << (i0+1) << " "  << (1+Th.Gh.Number(Th.edges[i0].on)) <<  endl;
    900                 if (Th.NbCrackedEdges)
    901         {
    902           f << "\nCrackedEdges\n"<< Th.NbCrackedEdges << endl;   
    903           for(i0=0;i0< Th.NbCrackedEdges; i0++)
    904             {
    905               f << Th.Number(Th.CrackedEdges[i0].a.edge) << " " ;
    906               f  << Th.Number(Th.CrackedEdges[i0].b.edge) << endl;
    907             }
    908         }
    909    } 
    910    if (&Th.BTh != &Th && Th.BTh.OnDisk && Th.BTh.name)
    911      {
    912        int *mark=new int[Th.nbv];
    913        Int4 i;
    914        for (i=0;i<Th.nbv;i++)
    915          mark[i]=-1;
    916        f << "\nMeshSupportOfVertices\n" <<endl;
    917        WriteStr(f,Th.BTh.name);
    918        f <<endl;
    919        f << "\nIdentityOfMeshSupport" << endl;
    920        WriteStr(f,Th.BTh.identity);
    921        f<<endl;
    922 
    923        f << "\nVertexOnSupportVertex" << endl;
    924        f<< Th.NbVertexOnBThVertex << endl;
    925        for(i=0;i<Th.NbVertexOnBThVertex;i++) {
    926          const VertexOnVertex & vov = Th.VertexOnBThVertex[i];
    927          Int4 iv = Th.Number(vov.v);
    928          mark[iv] =0;
    929          f << iv+1<< " " << Th.BTh.Number(vov.bv)+1 << endl;}
    930 
    931        f << "\nVertexOnSupportEdge" << endl;
    932        f << Th.NbVertexOnBThEdge << endl;
    933        for(i=0;i<Th.NbVertexOnBThEdge;i++) {
    934          const VertexOnEdge & voe = Th.VertexOnBThEdge[i];
    935          Int4 iv = Th.Number(voe.v);
    936          //      assert(mark[iv] == -1]);
    937          mark[iv] = 1;
    938          f << iv+1 << " " << Th.BTh.Number(voe.be)+1 << " " << voe.abcisse <<  endl;}
    939        
    940        f << "\nVertexOnSupportTriangle" << endl;   
    941        Int4 k = Th.nbv -  Th.NbVertexOnBThEdge - Th.NbVertexOnBThVertex;
    942        f << k << endl;
    943        //       Int4 kkk=0;
    944        CurrentTh=&Th.BTh;
    945        for (i=0;i<Th.nbv;i++)
    946          if (mark[i] == -1) {
    947            k--;
    948            Icoor2 dete[3];
    949            I2 I = Th.BTh.toI2(Th.vertices[i].r);
    950            Triangle * tb = Th.BTh.FindTriangleContening(I,dete);
    951            if (tb->link) // a true triangle
    952              {
    953                Real8 aa= (Real8) dete[1]/ tb->det, bb= (Real8) dete[2] / tb->det;
    954                f << i+1 << " " << Th.BTh.Number(tb)+1 << " " << aa << " " << bb << endl ;
    955              }
    956            else
    957              {
    958                double aa,bb,det[3];
    959                TriangleAdjacent ta=CloseBoundaryEdgeV2(I,tb,aa,bb);
    960                int k = ta;
    961                det[VerticesOfTriangularEdge[k][1]] =aa;
    962                det[VerticesOfTriangularEdge[k][0]] = bb;
    963                det[OppositeVertex[k]] = 1- aa -bb;
    964                Triangle * tb = ta;
    965                f << i+1 << Th.BTh.Number(tb)+1 << " " << det[1] << " " << det[2] <<endl;
    966              }
    967          }
    968        assert(!k);
    969        delete [] mark;
    970          
    971 
    972      }
    973    f << "\nEnd" << endl;
    974    //  Th.ConsLinkTriangle();
    975    delete [] reft;
    976    return f;
    977    
    978 }
    979 
    980 
    981216
    982217void Geometry::Write(const char * filename)
Note: See TracChangeset for help on using the changeset viewer.