bamgfreefem.cpp

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

      t=a;
      double delta = (b-a)/n;
      for ( nn=0;nn<=n;nn++,i++, t += delta)
        {
          if (nn==n) t=b; // to remove roundoff error 
          mp.label = k->label();
          k->code(stack); // compute x,y, label
          // cout << " ----- " << i << " " << mp.P.x << " " << mp.P.y << endl;
          vertices[i].r.x=mp.P.x;
          vertices[i].r.y=mp.P.y;
          vertices[i].ReferenceNumber=  mp.label;
          vertices[i].color = i;
          if (nn>0) {
            lmin=min(lmin,Norme2_2( vertices[i].r-vertices[i-1].r));
          }     
        }
    }
  lmin = sqrt(lmin);
  double eps = (lmin)/16.; 
  int nbvprev = i;
  long nbv=0;
  Gh->pmin =  vertices[0].r;
  Gh->pmax =  vertices[0].r;
  // recherche des extrema des vertices pmin,pmax
  for (i=0;i<nbvprev;i++) 
    {
      Gh->pmin.x = Min(Gh->pmin.x,vertices[i].r.x);
      Gh->pmin.y = Min(Gh->pmin.y,vertices[i].r.y);
      Gh->pmax.x = Max(Gh->pmax.x,vertices[i].r.x);
      Gh->pmax.y = Max(Gh->pmax.y,vertices[i].r.y);
    }
    
  double diameter=Max(Gh->pmax.x-Gh->pmin.x,Gh->pmax.y-Gh->pmin.y);
  Gh->coefIcoor= (MaxICoor)/diameter;
  Icoor1 epsI = (Icoor1) (Gh->coefIcoor*eps);
  ffassert(Gh->coefIcoor >0);

 if(lmin<diameter*1e-7) {
    ExecError(" Error points  border points to close < diameter*1e-7 ");}
  
  if (verbosity>2) 
    {
      cout <<"\t\t"  << "     Geom: min="<< Gh->pmin << "max ="<< Gh->pmax 
           << " hmin = " << Gh->MinimalHmin() <<  endl;
    }
  nbv = 0;        
  {  // find common point 
    QuadTree quadtree;
    Metric Id(1.);
    for ( i=0;i<nbvprev;i++) 
      { 
        vertices[i].i = Gh->toI2(vertices[i].r);
        vertices[i].m = Id;
        Vertex *v= quadtree.ToClose(vertices[i],eps,epsI,epsI) ;
        // quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 
        if( v && Norme1(v->r - vertices[i]) < eps )
          { vertices[i].color=v->color; }
        else  {quadtree.Add(vertices[i]);
        vertices[i].color = nbv++;}
      }
    /*   
         if (nbvprev-nbv==0)
         {
         reffecran();
         bamg::R2 O((Gh->pmin+Gh->pmax)/2),D(Gh->pmax-Gh->pmin);
         cadreortho(O.x,O.y,Max(D.x,D.y)*1.1);
         xGrafCoef = Gh->coefIcoor;
         yGrafCoef = Gh->coefIcoor;
         xGrafOffSet = Gh->pmin.x;
         yGrafOffSet = Gh->pmin.y;  
         quadtree.Draw();
         for (int i=0;i<nbvprev;i++)
         {
         rmoveto(vertices[i].r.x,vertices[i].r.y);
         
         char buf[100];
         sprintf(buf,"%d",i);
         plotstring(buf);
         }
         rattente(1);
         }
    */   
  } // to delete quadtree
  if (verbosity>1) 
  cout << " Nb of common points " << nbvprev-nbv <<endl;
  
  Gh->nbvx = nbv;
  Gh->nbv = nbv;
  
  Gh->vertices = new GeometricalVertex[nbv];
  throwassert(Gh->nbvx >= Gh->nbv);
  Gh->nbiv = Gh->nbv;
  // Int4 k=0;
  const Direction NoDirOfSearch;
  //  compression of points    
  int kkk;     
  for ( i=0,kkk=0;kkk<nbvprev;kkk++) 
    {
      if (vertices[kkk].color == i) // if new points 
        {
          Gh->vertices[i].r.x = vertices[kkk].r.x ;
          Gh->vertices[i].r.y = vertices[kkk].r.y;
          //Gh->vertices[i].link = Gh->vertices + i;
          throwassert(Gh->vertices[i].IsThe());
          Gh->vertices[i].ReferenceNumber = vertices[kkk].ReferenceNumber  ;
          Gh->vertices[i].DirOfSearch = NoDirOfSearch;
          Gh->vertices[i].color =0;
          Gh->vertices[i].Set();
          //  vertices[i].SetCorner();
	    if(Requiredboundary)
           Gh->vertices[i].SetRequired();
          i++;
        }
    }
  throwassert(i==nbv);
  R2 zero2(0,0);
  if(verbosity>5) 
    cout <<"\t\t"  << "     Record Edges: Nb of Edge " << Gh->nbe <<endl;
  throwassert(Gh->edges);
  throwassert (Gh->nbv >0); 
  Real4 *len =0;
  if (!hvertices) 
    {
      len = new Real4[Gh->nbv];
      for(i=0;i<Gh->nbv;i++)
        len[i]=0;
    }
  int nnn=0;
  i=0;
  for (E_BorderN const * k=b;k;k=k->next)
    {
      double & t = *  k->var(stack);
      double a(k->from(stack)),b(k->to(stack));
      n=Max(Abs(k->Nbseg(stack)),1L);     
      double delta = (b-a)/n;
      t=a+delta/2; 
      for ( nn=0;nn<n;nn++,i++, t += delta)
        {
          
          mp.label = k->label();
          k->code(stack); 
          Int4 i1 =  vertices[nnn].color, i2 =  vertices[++nnn].color;
          throwassert(i1 >= 0 && i1 < nbv);
          throwassert(i2 >= 0 && i2 < nbv);
          Gh->edges[i].ref = mp.label;
          Gh->edges[i].v[0]=  Gh->vertices + i1;
          Gh->edges[i].v[1]=  Gh->vertices + i2;
          R2 x12 = Gh->vertices[i2].r-Gh->vertices[i1].r;
          Real8 l12=Norme2(x12);
          Gh->edges[i].tg[0]=zero2;
          Gh->edges[i].tg[1]=zero2;
          Gh->edges[i].SensAdj[0] = Gh->edges[i].SensAdj[1] = -1;
          Gh->edges[i].Adj[0] = Gh->edges[i].Adj[1] = 0;
          Gh->edges[i].flag = 0;
          Gh->edges[i].link=0;
	  if(Requiredboundary)
	  Gh->edges[i].SetRequired();
	    
          if (!hvertices) 
            {
              Gh->vertices[i1].color++;
              Gh->vertices[i2].color++;
              len[i1] += l12;
              len[i2] += l12;
            }
          
          Hmin = Min(Hmin,l12);
        }
      nnn++; 
    }
  
  delete [] vertices; vertices=0;
  
  throwassert(nnn==nbvprev);
  throwassert(i==Gh->nbe);
  // definition  the default of the given mesh size 
  if (!hvertices) 
    {
      for (i=0;i<Gh->nbv;i++) 
        if (Gh->vertices[i].color > 0) 
          Gh->vertices[i].m=  Metric(len[i] /(Real4) Gh->vertices[i].color);
        else 
          Gh->vertices[i].m=  Metric(Hmin);
      delete [] len;
      
      if(verbosity>3) 
        cout <<"\t\t"  << "     Geom Hmin " << Hmin << endl;
    }
  
  Gh->NbSubDomains=nbsd;
  if (Gh->NbSubDomains>0) 
    {
      Gh->subdomains = new GeometricalSubDomain[  Gh->NbSubDomains];
      Int4 i1=0;
      i=0;
      for (E_BorderN const * k=b;k;k=k->next,i++)
        { long Nbseg =k->Nbseg(stack);
        long n=  Max(1L,Abs(Nbseg)); 
        Gh->subdomains[i].sens = Nbseg >0 ? 1 : -1;
        Gh->subdomains[i].edge=Gh->edges + i1;
        Gh->subdomains[i].ref = i;
        i1 += n;
        }
    }
  Gh->NbEquiEdges=0;
  Gh->NbCrackedEdges=0;
  Fem2D::Mesh * m=0;
  if (justboundary)
    m=bamg2msh(*Gh);
  else 
  {
      Gh->AfterRead();  
      int nbtx= nbvmax ? nbvmax :  (Gh->nbv*Gh->nbv)/9 +1000;
      
      Triangles *Th = 0;
      try { 
	  Th =new Triangles( nbtx ,*Gh);
      }
      catch(...)
      {
	  Gh->NbRef=0; 
	  delete Gh;
	  cout << " catch Err bamg "  << endl;
	  throw ;
     }
      m=bamg2msh(Th,true);      
      delete Th;
  }
  
  delete Gh;
  /* deja fait  dans bamg2msh
     Fem2D::R2 Pn,Px;
     m->BoundingBox(Pn,Px);
     m->quadtree=new Fem2D::FQuadTree(m,Pn,Px,m->nv);
  ---------- */
  mp=mps;   
 // m->decrement();
 //   Add2StackOfPtr2FreeRC(stack,m);// fait au niveau d'apres  07/2008 FH 
    
   
  return m;   
}

void E_BorderN::BoundingBox(Stack stack,double  &xmin,double & xmax, double & ymin,double & ymax) const
{  
  Fem2D::MeshPoint & mp (*Fem2D::MeshPointStack(stack)), mps = mp;
  for (E_BorderN const * k=this;k;k=k->next)
    {
      assert(k->b->xfrom); // a faire 
      double & t = *  k->var(stack);
      double a(k->from(stack)),b(k->to(stack));
      long n=Max(Abs(k->Nbseg(stack)),1L);     
      t=a;
      double delta = (b-a)/n;
      for (int  nn=0;nn<=n;nn++, t += delta)
        {
          if (nn==n) t=b; // to remove roundoff error 
          mp.label = k->label();
          k->code(stack); // compute x,y, label
          xmin=Min(xmin,mp.P.x);
          xmax=Max(xmax,mp.P.x);
          ymin=Min(ymin,mp.P.y);
          ymax=Max(ymax,mp.P.y);
        }
    }
  mp=mps; 
}  
void E_BorderN::Plot(Stack stack) const
{  
  using Fem2D::R2;
  
  Fem2D::MeshPoint & mp (*Fem2D::MeshPointStack(stack)), mps = mp;
  float x0,x1,y0,y1;
  getcadre(x0,x1,y0,y1);
  float h= (x1-x0)*0.01;
  int nbd=0;
  for (E_BorderN const * k=this;k;k=k->next)
    {
      nbd++;
      assert(k->b->xfrom); // a faire 
      double & t = *  k->var(stack);
      double a(k->from(stack)),b(k->to(stack));
      long n=Max(Abs(k->Nbseg(stack)),1L);     
      t=a;
      double delta = (b-a)/n;
      R2 P,Po;
      for (int  nn=0;nn<=n;nn++, t += delta)
        {
          if (nn==n) t=b; // to remove roundoff error 
          mp.label = k->label();
          mp.P.z=0;
          k->code(stack); // compute x,y, label
          P=mp.P.p2();
          couleur(2+mp.label);
          if(nn!=0) { LineTo(P);
          R2 uv(Po,P);
          double l = Max(sqrt((uv,uv)),1e-20);
          
          R2 dd = uv*(-h/l);
          R2 dn = dd.perp()*0.5;
          
          LineTo(P+dd+dn);
          MoveTo(P+dd-dn);
          LineTo(P);}
          else {
            DrawMark(mp.P.p2(),0.01);
            MoveTo(mp.P.p2());
            
          }
          
          
          //  cout << k->label()<< " " << nn << ", x,y = " << mp.P.x  << " , " << mp.P.y << endl;
          Po=P;
        }
      DrawMark(mp.P.p2(),0.01);
      MoveTo(mp.P.p2());
    }
  if(verbosity>9) cout << " -- Plot size : " << nbd << " Border \n";
  mp=mps; 
}
void E_BorderN::SavePlot(Stack stack,PlotStream & plot ) const
{  
    using Fem2D::R2;
    
    Fem2D::MeshPoint & mp (*Fem2D::MeshPointStack(stack)), mps = mp;
    //float x0,x1,y0,y1;
    //getcadre(x0,x1,y0,y1);
    //float h= (x1-x0)*0.01;
    
    long nbd1=0;
    for (E_BorderN const * k=this;k;k=k->next)
	nbd1++;
   plot << nbd1;
   int nbd=0;
    for (E_BorderN const * k=this;k;k=k->next)
      {
	  nbd++;
	  assert(k->b->xfrom); // a faire 
	  double & t = *  k->var(stack);
	  double a(k->from(stack)),b(k->to(stack));
	  long n=Max(Abs(k->Nbseg(stack)),1L);     
	  t=a;
	  double delta = (b-a)/n;
	  R2 P,Po;
	  plot<< (long) n;
	  for (int  nn=0;nn<=n;nn++, t += delta)
	    {
		if (nn==n) t=b; // to remove roundoff error 
		mp.label = k->label();
		mp.P.z=0;
		k->code(stack); // compute x,y, label
		P=mp.P.p2();
		plot << (long) mp.label <<P.x << P.y;
	    }
	 
      }
    assert(nbd==nbd1);
    if(verbosity>9) cout << " -- Plot size : " << nbd << " Border \n";
    mp=mps; 
}

Fem2D::Mesh *  BuildMeshBorder(Stack stack, E_BorderN const * const & b) 
{
  return BuildMesh(stack,b,true,0,true);
}
Fem2D::Mesh *  BuildMesh(Stack stack, E_BorderN const * const & b,bool Requiredboundary) 
{
  return BuildMesh(stack,b,false,0,Requiredboundary);
}

Fem2D::Mesh *  ReadTriangulate( string  * const & s) {
  using namespace Fem2D;
  KN<R2> xy;
  char c;
  int nv;
  for(int step=0;step<2;step++)
    {
      nv=0;
      ifstream f(s->c_str());
      if(!f) {cerr <<" Error openning file " << *s << endl;
      ExecError("Openning file ");}
      while (f.good())
        {
          R2 P;
          f >> P ;
          if (!f.good()) break;
          if (step) xy[nv]=P;
          nv++;
          while (f.get(c) &&  (c!='\n' && c!='\r' ) ) (void) 0; // eat until control (new line
        } 
      if (!step && nv ) xy.init(nv); // alloc the array        
    }
  if(verbosity)
    cout << " we read  " << nv << " coordinates  xy "<< endl;
  
  Mesh * m=new Mesh(nv,xy); 
  m->MakeQuadTree();
//  m->decrement(); //  07/2008 FH  auto del ptr
 //  delete s;  modif mars 2006 auto del ptr
  return m;
  
}
Fem2D::Mesh *  Triangulate( const  KN_<double> & xx,const  KN_<double> & yy) 
{
    using namespace Fem2D;
    ffassert(xx.N()==yy.N());
    int nv=xx.N();
    KN<R2> xy(nv);
    for(int i=0;i<nv;++i)
       xy[i]= R2(xx[i],yy[i]);
    Mesh * m=new Mesh(nv,xy); 
    m->MakeQuadTree();
   // m->decrement();  07/2008 FH  auto del ptr
    //  delete s;  modif mars 2006 auto del ptr
    return m;
    
}
Fem2D::Mesh *  ReadMeshbamg( string * const & s) {
  using bamg::Triangles;
  Triangles * bTh= new Triangles(s->c_str());
  // bTh->inquire();
  Fem2D::Mesh * m=bamg2msh(bTh,false);// no renum
  delete bTh;
  // delete s; modif mars 2006 auto del ptr
 //  m->decrement();  07/2008 FH  auto del ptr
  return m;
}

Fem2D::Mesh *  buildmeshbamg( string * const & s, int nbvxin) {

  using bamg::Triangles;
  using bamg::Geometry;
  Geometry Gh(s->c_str());
  double	 hmin = Gh.MinimalHmin();
  double	 hmax = Gh.MaximalHmax();
  int nbvx = nbvxin ? nbvxin : ((Gh.nbv*Gh.nbv)/9 +1000); 
  Triangles * bTh=  new Triangles(nbvx,Gh);
  // bTh->inquire();
  Fem2D::Mesh * m=bamg2msh(bTh,false);// no renum
  delete bTh;
 //  delete s; modif mars 2006 auto del ptr
 // m->decrement();
  return m;
}