inout_popenbinaire.c

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

#include "medit.h"
#include "libmesh5.h"
#include "extern.h"
#include "string.h"
#ifdef WIN32
#include <fcntl.h>
#endif

#define WrdSiz 4


void getline_bin_float_vertex(int ddim, double *c, int *ref){
  int     i;
  float   ff;
  for(i=0; i<ddim; i++ ){
    fread( (unsigned char *)&(ff) ,WrdSiz, 1, stdin);
    c[i] = ff;
  }
  fread( (unsigned char *)&(*ref) ,WrdSiz, 1, stdin);   
  

}

void getline_bin_double_vertex(int ddim, double *c, int *ref){
  int     i;
  
  for(i=0; i<ddim; i++ ) 
    fread( (unsigned char *)&(c[i]) ,WrdSiz, 2, stdin);   

  fread( (unsigned char *)&(*ref) ,WrdSiz, 1, stdin);   
  
}

void getline_bin_elem(int ddim, int *v, int *ref){
  int     i;
  
  for(i=0; i<ddim; i++){
    fread( (unsigned char *)&(v[i]) ,WrdSiz, 1, stdin);
  }
  fread( (unsigned char *)&(*ref) ,WrdSiz, 1, stdin);
}


void getline_bin_edge(int *v0, int *v1, int *ref){
  fread( (unsigned char *)&(*v0) ,WrdSiz, 1, stdin);
  fread( (unsigned char *)&(*v1) ,WrdSiz, 1, stdin);
  fread( (unsigned char *)&(*ref) ,WrdSiz, 1, stdin);
}


void getline_bin_int_noref(int ddim, int *v){
  int     i;
   for(i=0; i<ddim; i++){
    fread( (unsigned char *)&(v[i]) ,WrdSiz, 1, stdin);
  }
}

void getline_bin_float_noref(int ddim, double *v){
  int     i;
  float  ff;
  for(i=0; i<ddim; i++){
    fread( (unsigned char *)&(ff) ,WrdSiz, 1, stdin);
    printf("value of ff %f\n",ff);
    v[i] = ff;
  }
}

void getline_bin_double_noref(int ddim, double *v){
  int     i;
  for(i=0; i<ddim; i++){
    fread( (unsigned char *)&(v[i]) ,WrdSiz, 2, stdin);
  }
}


/**********************************/
/*   function for loadsol_popen   */

void read_TypeSizeTyptab_bin(int *type,int *size,int *typtab){

  char    data[256];
  char   *tictac;
  int     i;
  int    tmptype;
  int    tmpsize;

  fread( (unsigned char *)&(tmptype) ,WrdSiz, 1, stdin);
  assert( tmptype <=GmfMaxTyp );

  tmpsize=0;
  for(i=0; i<tmptype; i++){
    fread( (unsigned char *)&(typtab[i]) ,WrdSiz, 1, stdin);
    //printf("typtab[%i]=%i\n",i,typtab[i]);
    tmpsize += typtab[i];  
  }
  
  *size = tmpsize;
  *type = tmptype;
}

int loadMesh_popen_bin(pMesh mesh) {
  pPoint      ppt;
  pEdge       pr;
  pTriangle   pt;
  pQuad       pq;
  pTetra      ptet;
  pHexa       ph;
  double      d,dp1,dp2,dp3,dn[3];
  float      *n,fp1,fp2,fp3;
  int         i,ia,ib,inm,ref,is,k,disc,nn,nt,nq;
  char       *ptr,data[256];
  /* Rajout popen*/
  char       *tictac;
  char       *natureread;
  int        loopdebug;
  int        vatn[2];
  int        tvatn[3];
 
  /* variable rajouter pour popen */
  int retcode=0;
  int cod;
  int KwdCod;
  int NulPos;
  int NextKwdPos;
  /* rajout */
  float ftab[3];
#ifdef WIN32
  _setmode(fileno(stdin),O_BINARY);     
#endif 
  
  /* parse keywords */
  mesh->np   = 0;
  mesh->nt   = 0;
  mesh->nq   = 0;
  mesh->ntet = 0;
  mesh->nhex = 0;
  mesh->nc   = 0;
  mesh->nr   = 0;
  mesh->na   = 0;
  mesh->nri  = 0;
  mesh->nre  = 0;
  mesh->nvn  = 0;
  mesh->ntg  = 0;
  mesh->ne = mesh->nt + mesh->nq + mesh->ntet + mesh->nhex;


  loopdebug = -1;
  
  // read code
  fread( (unsigned char *)&cod ,WrdSiz, 1, stdin);
  if(cod != 1 ){
    printf("error in reading the binary file .meshb with popen\n");
    exit(1);
  }
  fread( (unsigned char *)&(mesh->ver) ,WrdSiz, 1, stdin);
  fread( (unsigned char *)&KwdCod ,WrdSiz, 1, stdin);
  if(KwdCod != GmfDimension ){
    printf("error in reading the binary file .meshb with popen\n");
    exit(1);
  }
  fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
  fread( (unsigned char *)&(mesh->dim) ,WrdSiz, 1, stdin);
  /*control of the dimension*/
  if( (mesh->dim != 2) && (mesh->dim != 3) ){
    printf("the dimension is not correct");
  }
  printf("reading dimension %i \n",mesh->dim); 
  while( !feof(stdin) ){

    loopdebug=loopdebug+1;
   
    fread( (unsigned char *)&KwdCod ,WrdSiz, 1, stdin);

    /* determination of KwdCod */
    printf("reading KwdCod %i %i\n",KwdCod,loopdebug);

    switch (KwdCod){
    
    case GmfVertices :
      natureread="Vertices";
      printf("reading %s",natureread);
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->np) ,WrdSiz, 1, stdin);
      printf(": number of vertices %i\n",mesh->np);
   
      if ( !mesh->np ) {
	fprintf(stdout,"  ## No vertex\n");
 	retcode=-1;
	goto Lret; 
      }
      if ( ddebug ) printf("allocate %d points\n",mesh->np);
      mesh->point = (pPoint)M_calloc(mesh->np+1,sizeof(Point),"zaldy1.point");
      assert(mesh->point);
      
      for (k=1; k<=mesh->np; k++) {
	//if(0) printf("lecture point du maillage k=%i np=%i ver=%i \n",k,mesh->np,mesh->ver);
	ppt = &mesh->point[k];
	
	if(mesh->ver==GmfFloat)
	  getline_bin_float_vertex(mesh->dim, ppt->c, &ref);
	else
	  getline_bin_double_vertex(mesh->dim, ppt->c, &ref);

	ppt->ref = ref & 0x7fff;
	ppt->tag = M_UNUSED;
	
      }
      ppt = &mesh->point[k];

      break;
    case GmfTriangles :
      natureread = "Triangles";
      printf("reading %s",natureread);
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nt) ,WrdSiz, 1, stdin);
      printf(": number of triangles %i\n",mesh->nt);
      if ( ddebug ) printf("allocate %d tria\n",mesh->nt);
      mesh->tria = (pTriangle)M_calloc(mesh->nt+1,sizeof(Triangle),"zaldy1.tria");
      assert(mesh->tria);

      disc = 0;
      for (k=1; k<=mesh->nt; k++) {
	pt = &mesh->tria[k];
	getline_bin_elem( 3, pt->v, &ref);
	
	pt->ref  = ref & 0x7fff;
	for (i=0; i<3; i++) {    
	  if ( pt->v[i] < 1 || pt->v[i] > mesh->np ) {
	    disc++;
	    pt->v[0] = 0;
	    break;
	  }
	  else {
	    ppt = &mesh->point[pt->v[i]];
	    ppt->tag &= ~M_UNUSED;
	  }
	}
      }
      break;
    case GmfQuadrilaterals :
      natureread = "Quadrilateral";
      printf("reading %s",natureread);
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nq) ,WrdSiz, 1, stdin);
      printf(": number of hexahedrons %i\n",mesh->nq);

      if ( ddebug ) printf("allocate %d quad\n",mesh->nq);
      if ( mesh->nq ) {
	mesh->quad = (pQuad)M_calloc(mesh->nq+1,sizeof(Quad),"zaldy1.quad");
	assert(mesh->quad);
      }

      for (k=1; k<=mesh->nq; k++) {
	pq = &mesh->quad[k];
	getline_bin_elem( 4, pq->v, &ref);
	pq->ref  = ref & 0x7fff;

	for (i=0; i<4; i++) {    
	  if ( pq->v[i] < 1 || pq->v[i] > mesh->np ) {
	    disc++;
	    pq->v[0] = 0;
	    break;
	  }
	  else {
	    ppt = &mesh->point[pq->v[i]];
	    ppt->tag &= ~M_UNUSED;
	  }
	}
      }
      
      break;
    case GmfTetrahedra :
      natureread = "Tetrahedra";
      printf("reading %s",natureread);
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->ntet) ,WrdSiz, 1, stdin);
      printf(": number of tetrahedrons %i\n",mesh->ntet);
      if ( mesh->ntet ) {
	if ( ddebug ) printf("allocate %d tetra\n",mesh->ntet);
	mesh->tetra = (pTetra)M_calloc(mesh->ntet+1,sizeof(Tetra),"zaldy1.tetra");
	assert(mesh->tetra);
      }

       for (k=1; k<=mesh->ntet; k++) {
	ptet = &mesh->tetra[k];
	getline_bin_elem( 4, ptet->v, &ref);
	ptet->ref  = ref & 0x7fff;
	for (i=0; i<4; i++) {    
	  if ( ptet->v[i] < 1 || ptet->v[i] > mesh->np ) {
	    disc++;
	    ptet->v[0] = 0;
	    break;
	  }
	  else {
	    ppt = &mesh->point[ptet->v[i]];
	    ppt->tag &= ~M_UNUSED;
	  }
	}
      }

      break;
    case GmfHexahedra :
      natureread = "Hexahedra";
      printf("reading %s",natureread);
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nhex) ,WrdSiz, 1, stdin);
      printf(": number of hexahedrons %i\n",mesh->nhex);

      if ( ddebug ) printf("allocate %d hexa\n",mesh->nhex);
      mesh->hexa = (pHexa)M_calloc(mesh->nhex+1,sizeof(Hexa),"zaldy1.hexa");
      assert(mesh->hexa);
            
      for (k=1; k<=mesh->nhex; k++) {
	ph = &mesh->hexa[k];
	getline_bin_elem( 8, ph->v, &ref);
	
	ph->ref  = ref & 0x7fff;
	for (i=0; i<8; i++) {    
	  if ( ph->v[i] < 1 || ph->v[i] > mesh->np ) {
	    disc++;
	    ph->v[0] = 0;
	    break;
	  }
	  else {
	    ppt = &mesh->point[ph->v[i]];
	    ppt->tag &= ~M_UNUSED;
	  }
	}
      }
      break;
    case GmfCorners :
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nc) ,WrdSiz, 1, stdin);
       
      for (k=1; k<=mesh->nc; k++) {
	fread( (unsigned char *)&is ,WrdSiz, 1, stdin);
	//fgets(data,256,stdin);  tictac = strtok(data," \n");  is = atoi(tictac);
	if ( is < 1 || is > mesh->np )
	  disc++;
	else {
	  ppt = &mesh->point[is];
	  ppt->tag |= M_CORNER;
	  ppt->tag &= ~M_UNUSED;
	}
      }
      break;
    case GmfRequiredVertices :
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nr) ,WrdSiz, 1, stdin);
      for (k=1; k<=mesh->nr; k++) {
       	//fgets(data,256,stdin);  tictac = strtok(data," \n");  is = atoi(tictac);
	fread( (unsigned char *)&is ,WrdSiz, 1, stdin);
	if ( is < 1 || is > mesh->np )
	  disc++;
	else {
	  ppt = &mesh->point[is];
	  ppt->tag |= M_REQUIRED;
	  ppt->tag &= ~M_UNUSED;
	}
      }
      break;
    case GmfEdges :
      natureread="Edges";
      printf("reading %s",natureread);
      fread( (unsigned char *)&NulPos, WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->na), WrdSiz, 1, stdin);
      printf(": number of edges %i\n",mesh->na);
      if ( ddebug ) printf("allocate %d edges\n",mesh->na);
      mesh->edge = (pEdge)M_calloc(mesh->na+1,sizeof(Edge),"zaldy1.edge");
      assert(mesh->edge);

      for (k=1; k<=mesh->na; k++) {
	getline_bin_edge( &ia, &ib, &ref);
	if ( ia < 1 || ia > mesh->np || ib < 1 || ib > mesh->np )
	  disc++;
	else {
	  pr = &mesh->edge[k];
	  pr->v[0] = ia;
	  pr->v[1] = ib;
	  pr->ref  = ref & 0x7fff;
	  pr->tag  = !pr->ref ? M_NOTAG : M_TAG;
	  ppt = &mesh->point[ia];
	  ppt->tag &= ~M_UNUSED;
	  ppt = &mesh->point[ib];
	  ppt->tag &= ~M_UNUSED;
	}
      }
      break;
    case GmfRidges :
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nri) ,WrdSiz, 1, stdin);
      
      for (k=1; k<=mesh->nri; k++) {
	//getline_1int( natureread, &is);
	fread( (unsigned char *)&is, WrdSiz, 1, stdin);
	if ( is < 1 || is > mesh->na )
	  disc++;
	else {
	  pr = &mesh->edge[is];
	  pr->tag |= M_RIDGE;
	}
      }

      break;
    case GmfRequiredEdges :
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nre) ,WrdSiz, 1, stdin);

      for (k=2; k<=mesh->nre; k++) {
	//getline_1int( natureread, &is);
	fread( (unsigned char *)&is, WrdSiz, 1, stdin);
	if ( is < 1 || is > mesh->na )
	  disc++;
	else {
	  pr = &mesh->edge[is];
	  pr->tag |= M_REQUIRED;
	}
      }

      break;
    case GmfNormals :
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->nvn) ,WrdSiz, 1, stdin);

       /*  allocation  */ 
      if ( !mesh->ntg ) {
	mesh->extra = (pExtra)M_calloc(1,sizeof(Extra),"zaldy1.extra");
	assert(mesh->extra);
      }
      
      mesh->extra->n = (float*)M_calloc(3*mesh->nvn+1,sizeof(float),"inmesh");
      assert(mesh->extra->n);

      for (k=1; k<=mesh->nvn; k++) {
	n = &mesh->extra->n[3*(k-1)+1];	

	if(mesh->ver == GmfFloat)
	  getline_bin_float_noref( 3, dn);
	else
	  getline_bin_double_noref( 3, dn);

	n[0] = dn[0];
	n[1] = dn[1];
	n[2] = dn[2];
	
	d = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
	if ( d > 0.0 ) {
	  d = 1.0 / sqrt(d);
	  n[0] *= d;
	  n[1] *= d;
	  n[2] *= d;
	}
      }
    
      break;
    case GmfNormalAtVertices :
      if(!mesh->nvn){
	printf("The field Normal need to be reading before %s",natureread);
	exit(1);
      }
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->extra->iv) ,WrdSiz, 1, stdin);
      mesh->extra->nv = (int*)M_calloc(mesh->np+1,sizeof(int),"inmesh");
      assert(mesh->extra->nv);
      
      for (k=1; k<=mesh->extra->iv; k++) {
	//GmfGetLin(inm,GmfNormalAtVertices,&nn,&is);
	getline_bin_int_noref( 2, vatn);
	nn=vatn[0];
	is=vatn[1];
	if ( nn < 1 || nn > mesh->np )
	  disc++;
	else
	mesh->extra->nv[nn] = is;
      }
      break;
    case GmfNormalAtTriangleVertices :
      if(!mesh->nvn){
	printf("The field Normal need to be reading before %s",natureread);
	exit(1);
      }
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->extra->it) ,WrdSiz, 1, stdin);
      mesh->extra->nt = (int*)M_calloc(3*mesh->nt+1,sizeof(int),"inmesh");
      assert(mesh->extra->nt);
      
      for (k=1; k<=mesh->extra->it; k++) {
	//GmfGetLin(inm,GmfNormalAtTriangleVertices,&nt,&is,&nn);
	getline_bin_int_noref( 3, tvatn);
	tvatn[0] = nt;
	tvatn[1] = is;
	tvatn[2] = nn;