meshread.cpp

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

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//

/*
 
 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 ;
#endif

namespace 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() ;//