ffglut.cpp

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

	    if(what==6)
	      for(int sk=0;sk<nsubT;++sk)
		{
		  
		  for(int l=0;l<4;++l)
		    {
		      int iv= Ksub[lK++];
		      ii[l]= iv;
		      Pt[l]=Pn[iv];
		      ff[l]=v[o+iv];
		    }
		  
		  for(int i=0;i< plot.Viso.N();++i)
		    {
		      plot.color(i+4);			    
		      drawisoTet( Pt,ff,plot.Viso[i]);
		    }
		  
		}
	  }
        glEndList();  // fin de la list
      }
    win->unsetLighting();
    ShowGlerror("b mesh  OnePlotFE plot");  
    Plot(*Th,false,plot.drawmeshes,plot.drawborder,plot,gllists+2,&oklist[2]);
    ShowGlerror("OnePlotFE::Draw");
}    
void OnePlotFE::Draw(OneWindow *win)
{
  initlist();
  ThePlot & plot=*win->theplot;
  ShowGlerror("begin OnePlotFE plot");
  plot.SetDefIsoV();
  win->setLighting();
  //    OneWindow * win=plot.win;// bof bof  la struct est tres mauvaise . 
  assert(win);
  const Mesh & Th(*this->Th);
  int nsubT=NbOfSubTriangle(nsub);
  int nsubV=NbOfSubInternalVertices(nsub);
  int nK=v.N()/ Th.nt;
  if(debug>4)
  cout << "\t\t\tOnePlotMesh::Draw  " <<v.N() << " ,nt " << Th.nt << " " << nK << " " 
       << nsubV << " " << what << " ,nv " << Th.nv <<  endl;
  ffassert(v.N()== Th.nt*nK);
  ffassert(nK = nsubV*what);
  int o=0;
  KN<R2> Pn(nsubV);
  if((debug > 10)) cout << " " << nsubT << " " << nsubV << endl;
  
  if(plot.fill && what==1)
    glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
  else
    glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
  
  if(what==2)
    glDisable(GL_DEPTH_TEST);
  else 
    glEnable(GL_DEPTH_TEST);

  R coef = plot.coeff;
  double xmin,xmax,ymin,ymax;
  win->getcadre(xmin,xmax,ymin,ymax);
  double d= Max(ymax-ymin,xmax-xmin);
  R kk = 4*win->hpixel;
  R cc = win->hpixel*40;
  
  int klist=0;
  if( (what==1) ) 
    if ( plot.fill) klist=1; 
  //  cout << klist << " ... " << oklist[klist] << "  fill = " <<  plot.fill <<  endl;
  if (oklist[klist])
    glCallList(gllists+klist);
  else 
    {
      oklist[klist]=1;
      glNewList(gllists+klist,GL_COMPILE_AND_EXECUTE); // save  la list sans affichage
      for(int k=0;k<Th.nt;++k, o+= nK)
	{
	  const Mesh::Element & K=Th[k];
	  for(int i=0;i<nsubV;++i)
	    Pn[i]=K(SubInternalVertex(nsub,i));
	  if(what==1)
	    for(int sk=0;sk<nsubT;++sk)
	      {
		int i0=numSubTriangle(nsub,sk,0);
		int i1=numSubTriangle(nsub,sk,1);
		int i2=numSubTriangle(nsub,sk,2);
		
		R ff[3]={v[o+i0],v[o+i1],v[o+i2]};
		R2 Pt[3]={Pn[i0],Pn[i1],Pn[i2]};
		if(plot.fill)
		  plot.DrawIsoTfill( Pt, ff, plot.Viso,plot.Viso.N());
		else
		  plot.DrawIsoT( Pt, ff, plot.Viso,plot.Viso.N()); 
	      }
	  else // what ==2
	    for (int i=0,j=0;i<nsubV;++i)
	      {
		R2 P=Pn[i];
		R2 uv(v[o+j],v[o+j+1]);
		j+=2;
		R  l = Max(sqrt((uv,uv)),1e-30) ;
		int col = 2+dichotomie(plot.Varrow,l);
		if(debug>100) 
		  cout << uv << " l= " << l << " " << coef << " " <<col <<  endl;
		
		plot.color(2+col);
		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);	   
		glBegin(GL_LINES);          
		
		win->Seg(P,P+uv);
		
		if (l>kk) {
		  win->Seg(P+uv,P+uv+dd+dn);
		  win->Seg(P+uv,P+uv+dd-dn);
		}
		glEnd();		
	      }

	}
      glEndList();  // fin de la list
    }
  
  // if(plot.drawmeshes)
  //  if(what==2)
  //  glEnable(GL_DEPTH_TEST);  
  ShowGlerror("b mesh  OnePlotFE plot");  
  win->unsetLighting();
  Plot(Th,false,plot.drawmeshes,plot.drawborder,plot,gllists+2,&oklist[2]);
  ShowGlerror("OnePlotFE::Draw");
}

void OnePlotCurve::Draw(OneWindow *win)
{
  initlist();

  ThePlot & plot= *win->theplot;
  plot.SetColorTable(16) ;
    double z = plot.z0;
  
  glBegin(GL_LINE_STRIP);    
  plot.color(2);
  // cout << "nePlotCurve::Draw " << xx << " " << yy << endl;
  for (int i=0;i<xx.N();i++)
    {
      glVertex3d(xx[i],yy[i],z);
      
    }
  glEnd(); 
    
}

void OnePlotBorder::Draw(OneWindow *win)
{
  initlist();

  glDisable(GL_DEPTH_TEST);
  ThePlot & plot= *win->theplot;
  R h = 8*win->hpixel;
  
  double z = plot.z0;
  plot.SetColorTable(16) ; 
  
  // vector<vector<pair<long,R2> > > data;
  for(int i=0;i<data.size() ;++i)
    {
      vector<pair<long,R2> > & v=data[i];
      ShowGlerror("end OnePlotBorder::Draw  1");
      
      
      for(int j=1;j<v.size();++j)
	{
	  //	  cout <<v[j].first << endl;
	  plot.color(2+v[j].first);
	  R2 Po(v[j-1].second), Pn(v[j].second);
	  R2 uv(Po,Pn);
          double l = Max(sqrt((uv,uv)),1e-20);
          
          R2 dd = uv*(-h/l);
          R2 dn = dd.perp()*0.5;
	  glLineWidth(2); 
	  glBegin(GL_LINES);    
	  win->Seg(Po,Pn);
	  glEnd();
	  
	  glLineWidth(1);
	  glBegin(GL_LINES);    	  
	  if(j!=1)
	    {
	      win->Seg(Po,Po+dd+dn);
	      win->Seg(Po,Po+dd-dn);
	    }
	  glEnd();
	}
      
      ShowGlerror("end OnePlotBorder::Draw  2");
      
      glPointSize(7); 
      glBegin(GL_POINTS);    
      int l= v.size()-1;
      plot.color(2+v[0].first);
      glVertex3d(v[0].second.x,v[0].second.y,z);
      plot.color(2+v[l].first);
      glVertex3d(v[l].second.x,v[l].second.y,z);
      glEnd();
      glPointSize(1);
      ShowGlerror("end OnePlotBorder::Draw  3");
    }
  ShowGlerror("end OnePlotBorder::Draw");
  
}


OneWindow::OneWindow(int h,int w,ThePlot *p)
  :
  icurrentPlot(lplots.begin()), 
  lplotssize(0),
  height(h),width(w),theplot(0),hpixel(1),
  Bmin(0,0),Bmax(1,1),oBmin(Bmin),oBmax(Bmax),zmin(0),zmax(1),
  windowdump(false),help(false), rapz0(-1.),rapz(1),withlight(false)
{
  
  add(p);
}


void OneWindow::set(ThePlot *p)
{
  //ffassert(p);
    bool first = !theplot;
    bool change = theplot != p;
    theplot=p;
    if(p)
      {
	plotdim=p->plotdim;
      }
    rapz0 =-1; // to recompute the defalut rapz
    //    p->win=this;
    //    if(first)
    DefaultView() ;
    
    
}

void OneWindow::add(ThePlot *p)
{
  if(p) {
    lplots.push_back(p);
    lplotssize++;
    ++icurrentPlot;
    if(icurrentPlot==lplots.end())
      --icurrentPlot;// the previous
    if(icurrentPlot != lplots.end())
      set(*icurrentPlot);
    if( lplotssize>10)
      {
	bool isfirst = theplot == *lplots.begin();  
	cout << " delete a plot " << *lplots.begin() << endl;
	delete *lplots.begin();
	lplots.erase(lplots.begin());
        lplotssize--;
	if(isfirst) set(*lplots.begin()); // change to the next plot
      }
  }
  else 
    set(p);
}

void OneWindow::DefaultView() 
{
  if(theplot)
    {
      plotdim=theplot->plotdim;
      R3 A(theplot->Pmin),B(theplot->Pmax);
      R3 D(A,B);
      R dxy= max(D.x,D.y);
      zmax = theplot->fmax;
      zmin = theplot->fmin;
      theta=theplot->theta;
      phi=theplot->phi;
      if(theplot->datadim==3) rapz0=1;
      else   if(rapz0<=0)
	{ //  ( zmax-zmin )*rapz0 =  0.3 dxyy
	  rapz0  =  0.4* dxy/(zmax-zmin) ;
	  if(debug>0)
	    {
	      cout << " rapz0 = " << rapz0 ;
	      cout << " dz = " << zmax-zmin  << " dxy =" << dxy << endl;
	    }
	}
      rapz=rapz0;
      coef_dist=theplot->dcoef;
      focal=theplot->focal;      

      if(theplot->datadim==3)
	{
	  Bmin3=A;
	  Bmax3=B;
	}
      else
	{ // data plot 2d ou 1 d... 
	  if(theplot->boundingbox.size() ==4)
	    {
	      Bmin3.x=theplot->boundingbox[0];
	      Bmin3.y=theplot->boundingbox[1];
	      Bmax3.x=theplot->boundingbox[2];
	      Bmax3.y=theplot->boundingbox[3];	    
	    }
	  else 
	    {
	      Bmin3.x=A.x;
	      Bmin3.y=A.y;
	      Bmax3.x=B.x;
	      Bmax3.y=B.y;
	    }
	  Bmin3.z=theplot->fmin;
	  Bmax3.z=theplot->fmax;
	}
      Pvue3=(Bmin3+Bmax3)/2;
      
      
      
      D *=0.05;      
      if(theplot->boundingbox.size() !=4)
	{
	  A -= D;
	  B += D;
	}
      else
	{
	  R x1=theplot->boundingbox[0],y1=theplot->boundingbox[1];
	  R x2=theplot->boundingbox[2],y2=theplot->boundingbox[3];
	  A = R2(min(x1,x2),min(y1,y2));
	  B = R2(max(x1,x2),max(y1,y2));
	}
      
      if (theplot->aspectratio)
	cadreortho(A.p2(),B.p2());
      else 
	cadre(A.p2(),B.p2());
    }
  hpixel = (Bmax.x-Bmin.x)/width;
  
  // SetView() ;
}

void  OneWindow::SetScreenView() const
{

  glDisable(GL_TEXTURE_2D);
  glDisable(GL_DEPTH_TEST);
  glMatrixMode(GL_PROJECTION); 
  glLoadIdentity(); 
  glOrtho(0,width,0,height,-1,1);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

void  OneWindow::SetView()
{
  if(plotdim==3 && theplot)
    {
      glViewport(0, 0,width, height);
      
      glMatrixMode(GL_PROJECTION); 
      glLoadIdentity(); 
      R ratio= (double) width / (double)  height; 
      glMatrixMode(GL_MODELVIEW);
      glLoadIdentity(); 
      
      
      R aspect=ratio;
      R3 DD(Bmin3,Bmax3);
      DD.z *= rapz;
      R dmax= DD.norme();;
      R dist = 0.5*dmax/sin(focal/2)*coef_dist;
       cam.x=Pvue3.x+cos(phi)*cos(theta)*dist;
       cam.y=Pvue3.y+cos(phi)*sin(theta)*dist;
       cam.z=Pvue3.z*rapz+dist*sin(phi);  
      R znear=max(dist-dmax,1e-30);
      R zfare=dist+dmax;
      gluPerspective(focal*180./M_PI,aspect,znear,zfare);
      /*      
      if (eye)
	{
	  R dmm = -dmax*ceyes;
	  R dx = -dmm*sin(theta);
	  R dy = dmm*cos(theta);
	  camx += dx*eye;
	  camy += dy*eye;
	  }
      */
      if(debug>2)
	{
	  cout <<" setview 3d: rapz " <<  rapz << " cam: ";
	  cout << cam << " Pvue:" ;
	  cout << Pvue3  << " theta " << theta << "  phi = "<<  phi << endl;
	}
      gluLookAt(cam.x,cam.y,cam.z,Pvue3.x,Pvue3.y,Pvue3.z*rapz,0.,0.,1.);
      glScaled(1.,1.,rapz);   
      
    }
  else
    {
      ShowGlerror("Begin SetView");   
      glDisable(GL_DEPTH_TEST);
      glViewport(0, 0,width, height);
      R dz0,dz1,zm=0;
      if(plotdim==3)
	{
	    dz0=Bmin3.z;
	    dz1=Bmax3.z;
	}
      else
	{  
	    R zzmin = Min(zmin,theplot->fminT);
	    R zzmax = Max(zmax,theplot->fmaxT);    
	    R dz = (zzmax-zzmin);
	    zm=(zzmin+zzmax)*0.5;
	    //  to be sur  the the z size is no zero . 
	    dz = max(dz,(Bmax.x-Bmin.x)*0.1);
	    dz = max(dz,(Bmax.y-Bmin.y)*0.1);
	    dz0=-dz;
	    dz1 = dz;
	}
	
      if((debug>3 )) cout << "\t\t\t   SetView " << this << " " << Bmin  << " " 
			  << Bmax  << " dz  " << dz0 << " " << dz1  
			  << " theta " << theta << "  phi = "<<  phi << endl;
      ShowGlerror("0 Set MV");
      glMatrixMode(GL_MODELVIEW);
      glLoadIdentity();
      ShowGlerror(" Set MV");
      glMatrixMode(GL_PROJECTION); 
      glLoadIdentity(); 
      ShowGlerror(" Set PM 1");
      glOrtho(Bmin.x,Bmax.x,Bmin.y,Bmax.y,dz0,dz1);
      
      ShowGlerror(" Set PM 2");
      
      R2 M=(Bmin+Bmax)/2.;
      glTranslated(0,0,-zm);
      
      //glLineWidth(1);
      //glColor3d(0.,0.,0.);
      glGetDoublev(GL_PROJECTION_MATRIX,projMatrix);
      ShowGlerror(" Get PM");
      glGetDoublev(GL_MODELVIEW_MATRIX,modelMatrix);
      ShowGlerror(" Get MV");
      glGetIntegerv(GL_VIEWPORT,viewport);
      ShowGlerror(" Get VP");
  
      
      ShowGlerror("End SetView ");
  }
    
}
void  OneWindow::resize(int w,int h)
{  double ww=width,hh=height;
    width=w;
    height=h;
    if (theplot->aspectratio)
      {	
	  cadreortho(oBmin,oBmax);
      }
    
}

void  OneWindow::zoom(int w,int h,R coef)
{
    GLdouble x=w,y=height-h,z=(zmin+zmax)/2.;
    GLdouble xx,yy,zz;
    

    GLint ok= gluUnProject( x,y,z,modelMatrix,projMatrix,viewport,&xx,&yy,&zz);
    ShowGlerror(" UnPro .. ");
    if(debug>2)
      cout << " ok " << ok << " " << x << " " << y << " " << z 
	   << " -> " << xx << " " << yy << " " << zz << endl;
    R2  oD(oBmin,oBmax);
    R2  D(Bmin,Bmax);
    R2 O(xx,yy);// oBmin.x+D.x*xx/width,oBmin.y+D.y*yy/height); 
    if((debug > 3)) cout<< " zoom : "  << this << " O " << O  
			<< " " << coef << " D = "<<  D<< "as "
			<< theplot->aspectratio <<  endl;
    oD *= 0.5*coef;
    R2 A = O - oD;
    R2 B = O + oD;
    if (theplot->aspectratio)
	cadreortho(A,B);
    else 
	cadre(A,B);
}
void OneWindow::MoveXView(R dx,R dy) 
{
  R3 D(Bmin3,Bmax3);
  R3 dd( dx*D.x*sin(theta),-dx*D.y*cos(theta), - dy*D.z);
  if(debug>2)
  cout << " MoveXView  "<< dx << " " << dy << " " << D << " mn: " << Bmin3 <<"  mx :" << Bmax3 << " d=" << dd << endl;
  Pvue3 += dd/50.;
  // cout << xm << " " << ym << " " << zm << endl;
}

void OneWindow::cadre(R2 A,R2 B)
{
      
    oBmin=Bmin=A;
    oBmax=Bmax=B;
    hpixel = (Bmax.x-Bmin.x)/width;    
    
    
}

void OneWindow::getcadre(double &xmin,double &xmax,double &ymin,double &ymax)
{
    xmin =  Bmin.x;
    xmax =  Bmax.x;
    ymin = Bmin.y;
    ymax = Bmax.y;
    
}
void OneWindow::Display()
{
  ffassert(this && theplot);
  SetScreenView() ;
  glColor3d(0.,0.,0.);
     
  if(help)
    {
      theplot->DrawHelp(this);
      help=false;
    }
  else
    {
  ShowGlerror("Begin Display");
  
  //  SetView();
  if(theplot)
    theplot->Draw(this);
  ShowGlerror("After Display");
    }
}
void OneWindow::cadreortho(R2 A, R2 B)
{
    R2 D(A,B);
    oBmin=A;
    oBmax=B;
    
    double cxy =  D.y*width/ (D.x*height);
    
    if ( D.y*width < D.x*height)  
	// width -> infty => D.x la ref
	D.y = D.x*(double) height/ width;
    else // height -> infty => D.y la ref
	D.x = D.y*(double) width/height;
    R2 M=(A+B)/2., D2=D/2.;