ffglut.cpp

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

#ifdef __APPLE__
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif
//#include <pthread.h>
#include <limits>
#include <cfloat>
#include <cstdlib>
#include <cstdio>
#include <cmath>
using namespace std;

#include <fstream>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <vector>
#include <list>
#include <map>
#include <utility>

#include "rgraph.hpp"
#include "fem.hpp"
#include "RNM.hpp"
#include "Mesh3dn.hpp"

#include "PlotStream.hpp"

extern long verbosity;

// add for the gestion of the endianness of the file.
//PlotStream::fBytes PlotStream::zott; //0123;
//PlotStream::hBytes PlotStream::zottffss; //012345678;
// ---- FH

using namespace Fem2D;
using std::numeric_limits;
const R pi=M_PI;//4*atan(1.); 
using namespace std;

int debug=1;

#include "ffglut.hpp"

#include "ffthreads.hpp"



//Mutex MutexNextPlot;
Thread::Id tidRead=0;
bool NoMorePlot=false;
bool NoMorePlotTilte=false;
ThePlot *currentPlot=0, *nextPlot=0;
bool inThreadRead=false;
FILE *datafile=0;

static  bool TryNewPlot( void );


void LauchNextRead();
void WaitNextRead();
THREADFUNC(ThreadRead,fd);
//void * ThreadRead(void *fd);

int kread=-1;

int Fin(int code)
{
  WaitNextRead();
  if(!NoMorePlot && debug>2)
    cout << " exit before end  " << endl;
  exit(NoMorePlot ? 0  : 1);
}

int   ReadOnePlot(FILE *fp)
{ 
  int err=0;
   if(!fp) return -4; 
  err= feof(fp) ;
  if(err) return -2;
  err= ferror(fp) ;
  if(err) return -3;

  PlotStream f(fp);
  f.set_binary_mode();
  const char *  magic="#!ffglutdata2..";
  err=0;
  // init ..
  if(kread==-1)
    {
      for(int i=0;i<strlen(magic);i++)
	{ int c=getc(fp);
	  err += c != magic[i];
	  if(err) break;
	}
      if(err) {
	if(debug>2)
	 cout << " Err read heading " << endl;
	 goto Lreturn;
	 //return err;
	}
      kread++;
      if(debug>2) cout << " Read entete " << endl;
	  int c1 =getc(fp);//
      if(c1==13)	  
        int c2 =getc(fp);//	


    }
  long cas; 
  f >> cas;
  err=-1;
  if (feof(fp)) goto Lreturn ;
  if((debug > 2)) cout << " ReadOnePlot " << kread+1<< " cas = " << cas << " " << nextPlot << endl;
  if(cas==PlotStream::dt_newplot)  
    {
      assert(nextPlot==0);
      nextPlot = new ThePlot(f,currentPlot,++kread);
	if(debug>1)	
      cout << " next is build " << nextPlot<< " wait :" << nextPlot->wait << " -> " << kread <<  endl;
      assert(nextPlot);
      err=0;
    }
  else 
    {
      err=1;
      cout << " Error Cas inconnue (skip) " << endl;
    }
 Lreturn:
  f.set_text_mode();
  return err;
}


void TimerNextPlot(int value)
{
  // the routine to  until the end of nextplot.
  // we use gluttimerfunc functionnaly
  //  remark, if we miss we retry.
  // -----
  //  if(debug) cout << " TimeNextPlot  " << endl;
  value=min(1000,(value*3)/2);// try at leat every 1 second (not to heavy computation)
  if(TryNewPlot())
    glutPostRedisplay();
  else 
    glutTimerFunc(value,TimerNextPlot,value);
}

int SendForNextPlot()
{
  //  to send a event to plot the date sheet.
  // and out a timer to wait to the end of read..
  // every 25/ second..  = 1000/25 = 40 ms
  if(NoMorePlot)
    {
    if((debug > 0)) cout << " send signal For Next plot, skip: No More Plot !  " << endl;
    return 0;
    }
  if((debug > 1)) cout << " Try to read read  plot "<< endl;
  //  put a timer for wait to the  end of read
  glutTimerFunc(40,TimerNextPlot,40);
  
  return 1;
}

static  bool TryNewPlot( void )
{
  // the routine to try to see if the next plot is read or not. 
  // -----------------------------------------------------------
  bool ret=false;
  if(debug>2)
    cout << "  TryNewPlot   plot : " << currentPlot << " next = " << nextPlot << endl;;
  if (nextPlot!=0)
    {
      WaitNextRead();
      if(debug>1) cout << " change current plot to: " << nextPlot << " et  Lock Plot . " << endl;;

      AllWindows[glutGetWindow()]->add(nextPlot);
      //if(currentPlot) delete currentPlot; //  a change fait dans add 
      // MutexNextPlot.WAIT();      
      currentPlot=nextPlot;
      nextPlot=0;
      // MutexNextPlot.Free();
      LauchNextRead();  
      ret=true;
    }
  return ret;    
}



int  signep4(int i0,int i1,int i2,int i3)
{ // calcul du signe dans la permutation 
  int s =1;
  if(i0>i1) s=-s,Exchange(i0,i1);
  if(i1>i2) s=-s,Exchange(i1,i2);
  if(i2>i3) s=-s,Exchange(i2,i3); // i3 max
  if(i0>i1) s=-s,Exchange(i0,i1);
  if(i1>i2) s=-s,Exchange(i1,i2); // i2 max < i
  if(i0>i1) s=-s,Exchange(i0,i1);
  return s;
}
inline R3 bary(const R3 K[4],R f[4],int i0,int i1,R v)
{
  R d=f[i0]-f[i1];
  assert(fabs(d)>1e-20);
  R l1= (f[i0] - v)/ d;  //  == 1 si v = f[i1]  
  R l0 = 1. -l1;
  assert(l0 >=-1e-10 && l1 >= -1e-10);
  return K[i0]*l0 + K[i1]*l1; // == K[i1] si l1 ==1 => v = f[i1] 
}
void drawisoTet(const R3 K[4],R f[4],R v)
{
  static const int  nvfaceTet[4][3]  ={{3,2,1}, {0,2,3},{ 3,1,0},{ 0,1,2}}  ;//{ {2,1,3},{0,2,3},{1,0,3},{0,1,2} };

  R3 P[4];
  int nP=0;
  int np[4],nm[4];
  int km=0,kp=0;
  for (int i=0;i<4;++i)
    {
      if(f[i]<=v) nm[km++]=i;
      if(f[i]>=v) np[kp++]=i;
    }
  
  //cout << "km kp "<< km << " " << kp << endl;
  int h=-1,b[3];
  if(kp==1 && km==3)
    {
      h = np[0];
      b[0]=nvfaceTet[h][0];
      b[1]=nvfaceTet[h][1];
      b[2]=nvfaceTet[h][2];
    }
  if(km==1 && kp == 3)
    {
      h = nm[0];
      b[0]=nvfaceTet[h][0];
      b[2]=nvfaceTet[h][1];
      b[1]=nvfaceTet[h][2];
    }
  if(kp==2 && km==2)
    {//  cas quad 
      if(signep4(nm[0],nm[1],np[0],np[1]) < 0)
	Exchange(nm[0],nm[1]);
      //  le tet m[0],nm[1],np[0],np[1] est positif
      P[0]=bary(K,f,nm[0],np[0],v);
      P[1]=bary(K,f,nm[0],np[1],v);
      P[2]=bary(K,f,nm[1],np[1],v);
      P[3]=bary(K,f,nm[1],np[0],v);
      nP=4;      
    }
  else if (h>=0)
    { // cas triangle 
      P[0]=bary(K,f,h,b[0],v);
      P[1]=bary(K,f,h,b[1],v);
      P[2]=bary(K,f,h,b[2],v);
      nP=3;
    }
  

  /*
    if(nP)
    {
    cout << "+ " << np[0] << " - " << nm[0] << endl;
    cout << nP << " ;  ";
    for(int i=0;i<nP;++i)
    cout << P[i] << " ;  ";
    cout << endl;
    
    }
  */
    if(nP)
    {
      if(nP>2)
	{
	  R3 N(R3(P[0],P[1])^R3(P[0],P[2]));
	  N /= N.norme();
	  glNormal3d(N.x,N.y,N.z);
	}
      glBegin(GL_POLYGON);
      for(int i=0;i<nP;++i)
	glVertex3f(P[i].x, P[i].y,P[i].z); // 
      glEnd();
    }

  //  verification de l'orientation
  assert(nP < 3 || det(P[0],P[1],P[2],K[np[0]]) >=0)   ;
  assert(nP < 3 || det(P[0],P[1],P[2],K[nm[0]]) <=0)   ;
  
}


int dichotomie(const KN_<double>  &viso,R v) 
{
    int i=0,j=viso.N(),k;
    if  (v <viso[0] || v >viso[j-1]) 
	return -1;  
    while (i<j-1)    
	if ( viso[k=(i+j)/2]> v) j=k;
	else i=k;
    return i;
}


map<int,OneWindow *> AllWindows;
int  ShowGlerror(const char *s)
{
    GLint error = glGetError();
    if ( error != GL_NO_ERROR )
	printf("Attention %s erreur : %x \n",s,error);   
    return error;   
}


//  def des couleurs de la tables 
void DefColor(float & r, float & g, float & b,
              int k,int nb, bool hsv,bool grey,KN<R> colors)
{
    int nbcolors = colors.N()/3;
    if(k<=0) {  r=g=b=1.;} //  white
    else if (k==1)  { r=g=b=0.; } // black
    else if (k >= nb)   {  r=g=b=0.;} // black
    else if (grey) { float gg = 0.1+0.9*float(k-2)/(nb-3); r=g=b=gg;} 
    else if (nbcolors<=1) {  
	float h=float(k-2)/(nb-2),s=1.,v=1.;
	hsvToRgb(h,s,v,r,g,b); 
    return;}     
    else   { //  interpolation dans la table hsv    
	int i= (k-2); 
	int j0= i*(nbcolors-1) / (nb-2);
	int j1=j0+1;
	int i0=  j0*(nb-2)/(nbcolors-1);
	int i1=  j1*(nb-2)/(nbcolors-1);
	int j03=j0*3,j13=j1*3;
	float a=float(i1-i)/(i1-i0),a1=1-a;
	if (hsv)
	  {
	      float h = colors[j03+0]*a + colors[j13+0]*a1;
	      float s = colors[j03+1]*a + colors[j13+1]*a1;
	      float v = colors[j03+2]*a + colors[j13+2]*a1;
	  hsvToRgb(h,s,v,r,g,b); }
	else 
	  {
	      r = colors[j03+0]*a + colors[j13+0]*a1;
	      g = colors[j03+1]*a + colors[j13+1]*a1;
	      b = colors[j03+2]*a + colors[j13+2]*a1;
	  }
    }     
    
}

void Plot(const Mesh & Th,bool fill,bool plotmesh,bool plotborder,ThePlot & plot,GLint gllists,int * lok)
{
    glDisable(GL_DEPTH_TEST);

    ShowGlerror("begin Mesh plot");
    glPolygonMode(GL_FRONT_AND_BACK,GL_LINE); 
    R z1= plot.z0;
    R z2= plot.z0;
    
    
    double r=0,g=0,b=0;
    if((debug > 3)) cout<< " OnePlotMesh::Draw " << plotmesh << " " << plotborder << " " <<  Th.neb << " " << z1 << " "  << z2 << endl;
    // plot.SetColorTable(16) ; 
    bool cc[3]= { plotborder , plotmesh && fill , plotmesh };
    int kk=0;
    //for(int i=0;i<3;i++)
    //  cout << cc[i] << " " << lok[i] << " , ";
    //cout << endl;
    if(cc[kk])
      if(lok[kk])   glCallList(gllists+kk);
      else 
	{ 
	  lok[kk]=1;
	  glNewList(gllists+kk,GL_COMPILE_AND_EXECUTE ); // save  la list sans affichage
	  glLineWidth(2); 
	  glBegin(GL_LINES);    
	  for (int i=0;i<Th.neb;i++)
	    {
	      const BoundaryEdge & K(Th.be(i)); 
	      plot.color(1+abs(K.lab));
	      glVertex3d(K[0].x,K[0].y,z1);
	      glVertex3d(K[1].x,K[1].y,z1);
	      
	      
	    }
	  glEnd(); 	  
	  glLineWidth(1); 
	  glEndList();  // fin de la list
	}
      else ;
    
    kk++;	
    if(cc[kk])
      if(lok[kk])   glCallList(gllists+kk);
      else 
	{ 
	  lok[kk]=1;
	  glNewList(gllists+kk,GL_COMPILE_AND_EXECUTE ); // save  la list sans affichage
	  glPolygonMode(GL_FRONT,GL_FILL);//GL_FILL	
	  glBegin(GL_TRIANGLES);
	  for (int i=0;i<Th.nt;i++)
	    {
	      const Triangle & K(Th[i]);
	      plot.color(K.lab?1+abs(K.lab):0);
	      
	      //glColor3d(r,g,b);
	      int i0= Th(K[0]),  i1= Th(K[1]),   i2= Th(K[2]) ;    		
	      glVertex3d(K[0].x,K[0].y,z2);
	      glVertex3d(K[1].x,K[1].y,z2);
	      glVertex3d(K[2].x,K[2].y,z2);
	      
	    }    
	  glEnd();
	  glEndList();  //
	}
    
    kk++;
    if(cc[kk])
      if(lok[kk])   glCallList(gllists+kk);
      else 
	{ 
	  lok[kk]=1;
	  glNewList(gllists+kk,GL_COMPILE_AND_EXECUTE ); // save  la list sans affichage
	  glPolygonMode(GL_FRONT,GL_LINE);
	  glBegin(GL_TRIANGLES);
	  for (int i=0;i<Th.nt;i++)
	    {
		const Triangle & K(Th[i]);
		plot.color(fill? 1 : 1+abs(K.lab));
		int i0= Th(K[0]),  i1= Th(K[1]),   i2= Th(K[2]) ;    
		glVertex3d(K[0].x,K[0].y,z1);
		glVertex3d(K[1].x,K[1].y,z1);
		glVertex3d(K[2].x,K[2].y,z1);
		
	    }    
	  
	  glEnd();
	glEndList();  // fin de la list
      }

    ShowGlerror("end Mesh plot");

}


void Plot(const Mesh3 & Th,bool fill,bool plotmesh,bool plotborder,ThePlot & plot,GLint gllists,int * lok)
{
  typedef Mesh3::BorderElement BE;
  typedef Mesh3::Element Tet;
  glEnable(GL_DEPTH_TEST);
  /*
  if(fill)  glPolygonMode(GL_FRONT_AND_BACK,GL_LINE); 
  else glPolygonMode(GL_FRONT_AND_BACK,GL_LINE); 
  */
  ShowGlerror("begin Mesh plot");
  
  glPolygonMode(GL_FRONT_AND_BACK,GL_LINE); 
  
  R z1= plot.z0;
  R z2= plot.z0;
  
  double r=0,g=0,b=0;
  
  bool cc[3]= { plotborder , plotborder && fill , plotmesh };
  
  int kk=0;
  if(cc[kk])
    if(lok[kk])   glCallList(gllists+kk);
    else 
      { 
	lok[kk]=1;
	glNewList(gllists+kk,GL_COMPILE_AND_EXECUTE ); // save  la list sans affichage
	glLineWidth(1); 
	glBegin(GL_TRIANGLES);    
	for (int i=0;i<Th.nbe;i++)
	  {
	    const BE & K(Th.be(i)); 
	    plot.color(1+abs(K.lab));
	    R3 N(R3(K[0],K[1])^R3(K[0],K[2]));
	    N /= N.norme();
	    glNormal3d(N.x,N.y,N.z);
	    glVertex3d(K[0].x,K[0].y,K[0].z);
	    glVertex3d(K[1].x,K[1].y,K[1].z);
	    glVertex3d(K[2].x,K[2].y,K[2].z);
	  }
	glEnd(); 
	glLineWidth(1); 
	glEndList();  // fin de la list	  
      }
  
  kk++;
  ShowGlerror("end Mesh plot");
  
}



void OnePlotError::Draw(OneWindow *win)
{
  initlist();
  ThePlot & plot=*win->theplot;
  win->SetScreenView() ;
  glColor3d(0.,0.,0.);
  cout << " Error plot item empty " << item <<  endl;
  int i = 4;
  char s[100];
  sprintf(s,"Warning the item %d fot the plot is empty",item);
  win->Show(s,4+item*2);
  win->SetView() ;
}

void OnePlotMesh::Draw(OneWindow *win)
{
  initlist();
  ThePlot & plot=*win->theplot;
  Plot(*Th,plot.fill,true,true,plot,gllists,oklist);
  ShowGlerror("OnePlotMesh::Draw");
}
void OnePlotMesh3::Draw(OneWindow *win)
{
  initlist();
    ThePlot & plot=*win->theplot;
    Plot(*Th,plot.fill,true,true,plot,gllists,oklist);
    ShowGlerror("OnePlotMesh3::Draw");
}
void OnePlotFE3::Draw(OneWindow *win)
{
  initlist();
  
  ThePlot & plot=*win->theplot;
    ShowGlerror("begin OnePlotFE plot");
    plot.SetDefIsoV();
    if(plot.fill && what==6)
      glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
    else
      glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
    
    if(what==6)
      glEnable(GL_DEPTH_TEST);
    else 
      glEnable(GL_DEPTH_TEST);
    win->setLighting();
    if(oklist[0])
      glCallList(gllists+0);
    else
      { 
	oklist[0]=1;
        glNewList(gllists+0,GL_COMPILE_AND_EXECUTE); // save  la list sans affichage
	int nsubV=Psub.N();
	int nsubT=Ksub.N()/4;
	KN<R3> Pn(nsubV);
	int nK=v.N()/ Th->nt;
	for(int k=0,o=0;k<Th->nt;++k, o+= nK)
	  {
	    const Mesh3::Element & K=(*Th)[k];
	    int ii[4];// 
	    R ff[4];
	    R3 Pt[4];
	    for(int i=0;i<nsubV;++i)
	      Pn[i]=K(Psub[i]);
	    int lK=0;