meshread.cpp

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

     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