meshread.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

      
      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)
{
  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)
{

 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)
{
  //  
  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)
{
    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 )
{
  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;