drawing.cpp

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

  RN gk(nbdf);
  for (int i=0;i<nbdf;i++) // get the local value
    fk[i] = U[operator()(i)];
  for (int i=0;i<nbdf;i++) // get the local value
    gk[i] = V[operator()(i)];
    
  RNMK fb(nbdf,N,3); //  the value for basic fonction
  RN f0(3);
  RN f1(3);
  R2 Pt,P[3],A(T[0]),B(T[1]),C(T[2]);
  for (int k=0;k<nsb2;k++)
   {   
   //  ff=0.0;
     for (int j=0;j<3;j++)
      { 
     //   cout << nbdf << endl;
        Pt=SubTriangle(nsb,k,j); //  current point 
        BF(whatd,Pt,fb);
        f0[j] = (fb('.',i0,0),fk);
        f1[j] = (fb('.',i1,0),gk);
     //   if(number<2) 
     //   cout << number << " " << fk << " " << fb('.',0,0) << " :  " << f0[j] << " " <<  f1[j] << endl;
        P[j] = A*(1-Pt.x-Pt.y)+B*Pt.x + C*Pt.y;
        R2 uv(f0[j],f1[j]);
        R  l = Max(sqrt((uv,uv)),1e-30) ;
        
         {
          couleur(2+dichotomie(Viso,l));
          uv = coef*uv;
          l *= coef;
          R2 dd = uv*(-0.005/l);
          R2 dn = dd.perp()*0.5;
          if (l*10000.< kk) continue;
          if (l < kk) 
	        uv = uv*(kk/l);
          else if (l> cc)
	        uv = uv*(cc/l);	   
          MoveTo(P[j]);       
          LineTo(P[j]+uv);
          if (l>kk) {
            LineTo(P[j]+uv+dd+dn);
            MoveTo(P[j]+uv+dd-dn);
            LineTo(P[j]+uv);}
          }
      }
   } 
} 

void DrawIsoT(const R2 Pt[3],const R ff[3],const RN_ & Viso)
{
  R2 PQ[5];
  int NbIso = Viso.N();
  for(int l=0;l< NbIso;l++)  /*    loop on the level curves */
    {
      R xf = Viso[l];
      int im=0;
      for(int i=0;i<3;i++) // for the  3 edges 
	{
	  int j = (i+1)%3;
	  R fi=(ff[i]);
	  R fj=(ff[j]);
	  
	  if(((fi<=xf)&&(fj>=xf))||((fi>=xf)&&(fj<=xf)))
	    {
	      if (Abs(fi-fj)<=0.1e-10) 	/* one side must be drawn */
		{
		  couleur(l+4);
		  MoveTo(Pt[i]);
		  LineTo(Pt[j]);
		}
	      else
		{
		  R  xlam=(fi-xf)/(fi-fj);
		  PQ[im++]   = Pt[i] * (1.F-xlam)  +  Pt[j]* xlam;
		}
	    }
	}
      
      if (im>=2) /*    draw one segment */
	{
	  couleur(l+4);
	  MoveTo(PQ[0]);
	  LineTo(PQ[1]);
	}
    }
  
} 
void DrawIsoTfill(const R2 Pt[3],const R ff[3],const RN_ & Viso,double rapz)
{
  R2 PQ[10];
  int NbIso = Viso.N();
  R eps= (Viso[NbIso-1]-Viso[0])*1e-6;
  for(int l=1;l< NbIso;l++)  //   loop on the level curves 
    {
      R xfb = Viso[l-1];
      R xfh = Viso[l];
      int im=0;
      for(int i=0;i<3;i++) // for the  3 edges 
	{
          int j=(i+1)%3;
          R fi=(ff[i]);
          R fj=(ff[j]);
          R xxfb =  xfb;
          R xxfh =  xfh;
	  if (fj<fi ) Exchange(xxfb,xxfh);
          R xf  = xxfb;
	  if(((fi<=xf)&&(fj>=xf))||((fi>=xf)&&(fj<=xf)))
	    {
	      if (Abs(fi-fj)>=0.1e-20)
		{
		  R  xlam=(fi-xf)/(fi-fj);
		  PQ[im++]   = Pt[i] * (1.F-xlam)  +  Pt[j]* xlam;
		}
	    }
          xf = xxfh;	  
	  if(((fi<=xf)&&(fj>=xf))||((fi>=xf)&&(fj<=xf)))
	    {
	      if (Abs(fi-fj)>=0.1e-20)
		{
		  R  xlam=(fi-xf)/(fi-fj);
		  PQ[im++]   = Pt[i] * (1.F-xlam)  +  Pt[j]* xlam;
		}
	    }
	   if (  xfb-eps <=fj  && fj <= xfh+eps) 
	     PQ[im++] = Pt[j];
	}
       if (im>2) 
         {
           float f[20];
           int k=0;
           for (int i=0;i<im;i++)
            f[k++]=PQ[i].x,f[k++]=PQ[i].y;
           f[k++]=f[0];
           f[k++]=f[1];
           
           couleur(3+l);
           fillpoly(im,f);
         }
     }
} 

template<class R2>
void TBoundaryEdge<R2>::Draw() const 
{
  couleur(2+(lab%6));
  MoveTo(*vertices[0]);
  LineTo(*vertices[1]);
}

    void FElement::SaveDraw(const KN_<R> & U,int composante,R* Usave) const 
    {
	int nsb = nbsubdivision();
	int nsbv = NbOfSubInternalVertices(nsb);
	int nbdf=NbDoF();
	RN fk(nbdf);
	for (int i=0;i<nbdf;i++) // get the local value
	    fk[i] = U[operator()(i)];
	R2 Pt;
	bool whatd[last_operatortype];
	initwhatd(whatd,0);
	
	RNMK fb(nbdf,N,3); //  the value for basic fonction
	for (int k=0;k<nsbv;k++)
	  {
	      Pt=SubInternalVertex(nsb,k); 
	      BF(whatd,Pt,fb);
	      Usave[k] = (fb('.',composante,0),fk);
	  }
	
    }
    void FElement::SaveDraw(const KN_<R> & U,const KN_<R> & V,int iU,int iV,R * Usave) const 
    {
	int nsb = nbsubdivision();
	int nsbv = NbOfSubInternalVertices(nsb);
	int nbdf=NbDoF();
	RN fk(nbdf);
	RN gk(nbdf);
	for (int i=0;i<nbdf;i++) // get the local valu
	  {
	      fk[i] = U[operator()(i)];
	      gk[i] = V[operator()(i)];
	  }
	R2 Pt;
	bool whatd[last_operatortype];
	initwhatd(whatd,0);
	
	RNMK fb(nbdf,N,3); //  the value for basic fonction
	for (int k=0,k2=0;k<nsbv;k++)
	  {
	      Pt=SubInternalVertex(nsb,k); 
	      BF(whatd,Pt,fb);
	      Usave[k2++] = (fb('.',iU,0),fk);
	      Usave[k2++] = (fb('.',iV,0),gk);
	  }
	
    } 
    
KN<double>  FESpace::newSaveDraw(const KN_<R> & U,int composante,int & lg,int & nsb) const 
    {
	nsb = TFE[0]->nbsubdivision;
	int nsbv = NbOfSubInternalVertices(nsb);
	lg = nsbv*Th.nt;
	cout << "newSaveDraw what: nt " << Th.nt << " " << nsbv << " " << lg << endl;
	KN<double> v(lg);
	ffassert(v);
        for (int k=0,i=0;k<Th.nt;k++)
	  {
	    (*this)[k].SaveDraw( U,composante,&v[i]);	
	      i+=nsbv;
	  }
	return KN<double>(true,v);// to remove the copy.
    }
KN<double>  FESpace::newSaveDraw(const KN_<R> & U,const KN_<R> & V,int iU,int iV,int & lg,int & nsb) const 
    {
	nsb = TFE[0]->nbsubdivision;
	int nsbv = NbOfSubInternalVertices(nsb)*2;
	lg = nsbv*Th.nt;
	
	KN<double> v(lg);
	
        for (int k=0,i=0;k<Th.nt;k++)
	  {
	      (*this)[k].SaveDraw( U,V,iU,iV,&v[i]);	
	      i+=nsbv;
	  }
	return  KN<double>(true,v);// to remove the copy.
    }   
void  FESpace::Draw(const RN_& U,const RN_ & Viso,int j,float *colors,int nbcolors,bool hsv,bool drawborder) const 
{
  showgraphic();
  NewSetColorTable(Viso.N()+4,colors,nbcolors,hsv);
  for (int k=0;k<Th.nt;k++) 
    (*this)[k].Draw( U,Viso,j);
  NewSetColorTable(2+6,colors,nbcolors,hsv);
  if(drawborder) Th.DrawBoundary();
  NewSetColorTable(Viso.N()+4,colors,nbcolors,hsv);

}

void  FESpace::Drawfill(const RN_& U,const RN_ & Viso,int j,double rapz,float *colors,int nbcolors,bool hsv,bool drawborder) const 
{
  showgraphic();
  NewSetColorTable(Viso.N()+4,colors,nbcolors,hsv);
  for (int k=0;k<Th.nt;k++) 
    (*this)[k].Drawfill( U,Viso,j,rapz);
  NewSetColorTable(2+6,colors,nbcolors,hsv);
  if(drawborder) Th.DrawBoundary();
  NewSetColorTable(Viso.N()+4,colors,nbcolors,hsv);
}



R2 FESpace::MinMax(const KN_<R>& U,const KN_<R>& V,int j0,int j1,bool  bb) const 
{
  R2 Pminmax(1e100,-1e100);
    for (int k=0;k<Th.nt;k++)
     { 
        const Triangle & K(Th[k]);
        R2 A(K[0]),B(K[1]),C(K[2]);
        if (bb  || InRecScreen(Min(A.x,B.x,C.x),Min(A.y,B.y,C.y),Max(A.x,B.x,C.x),Max(A.y,B.y,C.y))) 
         Pminmax=minmax(Pminmax,(*this)[k].MinMax(U,V,j0,j1));
     }
   return Pminmax;
}
R2 FESpace::MinMax(const KN_<R>& U,int j0,bool bb) const 
{
  R2 Pminmax(1e100,-1e100);
    for (int k=0;k<Th.nt;k++) 
      {
        const Triangle & K(Th[k]);
        R2 A(K[0]),B(K[1]),C(K[2]);
        if (bb  || InRecScreen(Min(A.x,B.x,C.x),Min(A.y,B.y,C.y),Max(A.x,B.x,C.x),Max(A.y,B.y,C.y))) 
          Pminmax=minmax(Pminmax,(*this)[k].MinMax(U,j0));
      }
   return Pminmax;
}
void  FESpace::Draw(const KN_<R>& U,const RN_ & Viso, R coef,int j0,int j1,float *colors,int nbcolors,bool hsv,bool drawborder) const 
{ 
  showgraphic();
  NewSetColorTable(Viso.N()+5,colors,nbcolors,hsv);
  for (int k=0;k<Th.nt;k++) 
    (*this)[k].Draw( U,U,Viso,coef,j0,j1);
  NewSetColorTable(2+6,colors,nbcolors,hsv);
  if(drawborder) Th.DrawBoundary();
  NewSetColorTable(Viso.N()+5,colors,nbcolors,hsv);

}

void  FESpace::Draw(const KN_<R>& U,const KN_<R>& V,const RN_ & Viso, R coef,int iu,int iv,float *colors,int nbcolors,bool hsv,bool drawborder) const 
{ 
  showgraphic();
  NewSetColorTable(Viso.N()+5,colors,nbcolors,hsv);
  for (int k=0;k<Th.nt;k++) 
    (*this)[k].Draw( U,V,Viso,coef,iu,iv);
  NewSetColorTable(2+6,colors,nbcolors,hsv);
  if(drawborder) Th.DrawBoundary();
   NewSetColorTable(Viso.N()+5,colors,nbcolors,hsv);

}


template<class R2>
void TTriangle<R2>::Draw(double skrink) const
{
  const TTriangle<R2> & K(*this);
  R2 A(K[0]),B(K[1]),C(K[2]),G((A+B+C)/3.);
  A = G + (A-G)*skrink;
  B = G + (B-G)*skrink;
  C = G + (C-G)*skrink;
  MoveTo(A);
  LineTo(B);
  LineTo(C);
  LineTo(A);
}

template<class R2>
void TTriangle<R2>::Draw(int edge,double skrink) const 
{
  const TTriangle<R2> & K(*this);
  R2 A(K[0]),B(K[1]),C(K[2]),G((A+B+C)/3.);
  MoveTo(G+(*vertices[(edge+1)%3]-G)*skrink);
  LineTo(G+(*vertices[(edge+2)%3]-G)*skrink);
}

template<class R2>
void TTriangle<R2>::Fill(int color) const 
{
  const TTriangle<R2> & K(*this);
  R2 A(K[0]),B(K[1]),C(K[2]);
  float p[]={A.x,A.y,B.x,B.y,C.x,C.y};
  int c=LaCouleur(); 
  couleur(color); 
  fillpoly(3,p);
  couleur(c); 
}

void DrawMark(R2 P,R k)
 {
   float x0,x1,y0,y1;
   getcadre(x0,x1,y0,y1);
   float h= (x1-x0)*(float)k;
   rmoveto((float)P.x+h,(float)P.y-h);
   rlineto((float)P.x+h,(float)P.y+h);
   rlineto((float)P.x-h,(float)P.y+h);
   rlineto((float)P.x-h,(float)P.y-h);
   rlineto((float)P.x+h,(float)P.y-h);
 }
 
 Triangle * Mesh::Find(const R2 & P) const 
 {
    // brute force
    
    for (int i=0;i<nt;i++)
     {
      kthrough++;
      const Triangle & K(triangles[i]);
      R2 A(K[0]),B(K[1]),C(K[2]);
      R a=Area2(P,B,C);
      R b=Area2(A,P,C);
      R c=Area2(A,B,P);
      R s=a+b+c;
      R eps=s*1e-6;
      if (a>-eps && b >-eps && c >-eps) return triangles + i;
      
     }
   return 0; // outside 
 }

template<class R2>
 void TMortar<R2>::Draw() const {
   throwassert(Th);
   for (int i=0;i<nleft;i++)
     (*Th)[left[i]/3].Draw(left[i]%3,0.8);     
   for (int i=0;i<nright;i++)
     (*Th)[right[i]/3].Draw(right[i]%3,0.8);
     
 }
}