inout_popenbinaire.c

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

	if ( nt < 1 || nt > mesh->nt || is < 1 || is > 3 || nn < 1 || nn > mesh->nvn )
	  disc++;
	else
	  mesh->extra->nt[3*(nt-1)+is] = nn;
      }
      break;
    case GmfNormalAtQuadrilateralVertices :
      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->iq) ,WrdSiz, 1, stdin);
      mesh->extra->nq = (int*)M_calloc(4*mesh->nq+1,sizeof(int),"inmesh");
      assert(mesh->extra->nq);
      
      for (k=1; k<=mesh->extra->iq; k++) {
	//GmfGetLin(inm,GmfNormalAtQuadrilateralVertices,&nq,&is,&nn);
	getline_bin_int_noref( 3, tvatn);
	tvatn[0] = nq;
	tvatn[1] = is;
	tvatn[2] = nn;
	if ( nq < 1 || nq > mesh->nq || is < 1 || is > 4 || nn < 1 || nn > mesh->nvn )
	  disc++;
	else
	  mesh->extra->nq[3*(nq-1)+is] = nn;
      }
      break;
    case GmfTangents :
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->ntg) ,WrdSiz, 1, stdin);
      
      if ( !mesh->nvn ) {
	mesh->extra = (pExtra)M_calloc(1,sizeof(Extra),"zaldy1.extra");
	assert(mesh->extra);
      }
     
      mesh->extra->t = (float*)M_calloc(3*mesh->ntg+1,sizeof(float),"inmesh");
      assert(mesh->extra->t);

      for (k=1; k<=mesh->ntg; k++) {
	n = &mesh->extra->t[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 GmfTangentAtVertices :
      if( !mesh->ntg ){
	printf("The field Tangent need to be reading before %s",natureread);
	exit(1);
      }
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->extra->tv) ,WrdSiz, 1, stdin);

      mesh->extra->tv = (int*)M_calloc(mesh->np+1,sizeof(int),"inmesh");
      assert(mesh->extra->tv);

      for (k=1; k<=mesh->extra->jv; k++) {
	getline_bin_int_noref( 2, vatn);
	nn=vatn[0];
	is=vatn[1];
     
	if ( nn < 1 || nn > mesh->np )
	  disc++;
	else
	  mesh->extra->tv[nn] = is;
      }

      break;
    case GmfTangentAtEdgeVertices :
      if( !mesh->ntg ){
	printf("The field Tangent need to be reading before %s",natureread);
	exit(1);
      }
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&(mesh->extra->je) ,WrdSiz, 1, stdin);

      mesh->extra->te = (int*)M_calloc(2*mesh->na+1,sizeof(int),"inmesh");
      assert(mesh->extra->te);

      for (k=1; k<=mesh->extra->je; k++) {
	getline_bin_int_noref( 3, tvatn);
	nt=tvatn[0];
	is=tvatn[1];
	nn=tvatn[2];
	
	if ( nt < 1 || nt > mesh->np || is < 1 || is > 2 || nn < 1 || nn > mesh->ntg )
	  disc++;
	else
	  mesh->extra->te[3*(nt-1)+is] = is;
      }

      break;
    case GmfEnd :
     
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      printf("End of mesh\n");
      break;
    default :
      printf("This data KwdCod it is not taken in this version\n");
      exit(1);
      break;
    }
    
    if(KwdCod == GmfEnd){
      break;
    }
    /*
      if( KwdCod == GmfVertices ){
      }
      else if(KwdCod == GmfTriangles){
      }
      else if(KwdCod == GmfQuadrilaterals){
      }
      else if(KwdCod == GmfTetrahedra){
      }
      else if(KwdCod == GmfHexahedra){
      }
      else if(KwdCod == GmfCorners){
      }
      else if(KwdCod == GmfRequiredVertices){
      }
      else if(KwdCod == GmfEdges){
      }
      else if(KwdCod == GmfRidges){
      }
      else if(KwdCod == GmfRequiredEdges){
      }
      else if(KwdCod == GmfNormals){
      }
      else if(KwdCod == GmfTangents){
      }
      else if(KwdCod == GmfEnd){
      }
      else{
      printf("This data KwdCod it is not taken in this version\n");
      exit(1);
      }
    */
  }
  
  /* check if vertices and elements found */
  if ( !mesh->np ) {
    fprintf(stdout,"  ## No vertex\n");
    retcode=-1;
    goto Lret; 
  }
  mesh->ne = mesh->nt + mesh->nq + mesh->ntet + mesh->nhex;

  if ( disc > 0 ) {
    fprintf(stdout,"  ## %d entities discarded\n",disc);
  }
  retcode=1;

 Lret:
#ifdef WIN32
  _setmode(fileno(stdin),O_TEXT);     
#endif 
return retcode;
}

/* function of lecture */

int loadScaVecTen_bin(pMesh mesh, int numsol, int dim, int ver, int nel, int type, int size, int * typtab, int key) {
  pSolution    sol;
  double       dbuf[ GmfMaxTyp ];
  float        fbuf[ GmfMaxTyp ];
  double       m[6],lambda[3],eigv[3][3],vp[2][2];
  int          inm,k,i,iord,off;
  char        *ptr,data[128];
  double       ScaSol[1],VecSol[3],TenSol[9];
  float        fScaSol[1],fVecSol[3],fTenSol[9];
  int retcode=0;

  if(ddebug) printf("numsol=%i, type=%i, size=%i\n",numsol,type,size); 
  if ( numsol > type )  numsol = 1;
  numsol--;
  mesh->nbb    = nel;
  mesh->bbmin  =  1.0e20;
  mesh->bbmax  = -1.0e20;
  if ( !zaldy2(mesh) ) {
    retcode=0;
  }
  sol = mesh->sol;
  sol->dim = dim;
  sol->ver = ver;

  off = 0;
  for (i=0; i<numsol; i++) {
    if(ddebug) printf("typtab[%i]=%i",i,typtab[i]); 
    switch(typtab[i]) {
      case GmfSca:
        off++; break;
      case GmfVec:
        off += sol->dim; break;
      case GmfSymMat:
        off += sol->dim*(sol->dim+1)/2;  break;
    }
  }
  if(ddebug) printf("typtab[%i]=%i, off%i",i,typtab[i],off); 
  //printf("min= %f, max= %f\n",mesh->bbmin,mesh->bbmax);
 
  if(ddebug) printf("numsol=%i,typtab[i]=%i\n",numsol,typtab[i]); 
  fflush(stdout);
  switch(typtab[numsol]) {
  case GmfSca:
    mesh->nfield = 1;
    
    for (k=1; k<=nel; k++) {
      //   reading data must be a double !
      getline_bin_double_noref(1, ScaSol);
      mesh->sol[k].bb = ScaSol[0];
      
      if ( mesh->sol[k].bb < mesh->bbmin ) mesh->bbmin = mesh->sol[k].bb;
      if ( mesh->sol[k].bb > mesh->bbmax ) mesh->bbmax = mesh->sol[k].bb;
    }
    break;

  case GmfVec:
    if ( ddebug )  printf("   vector field\n");
    mesh->nfield = sol->dim;

    for (k=1; k<=nel; k++) {
      //   reading data must be a double ! 
      mesh->sol[k].bb = 0.0;      
      getline_bin_double_noref( sol->dim, VecSol);
      
      for (i=0; i<sol->dim; i++) {
	fbuf[off+i] = VecSol[i];
	//printf("solution vectorielle %i composante %i %f\n",k,i,VecSol[i]);
      }
      
      for (i=0; i<sol->dim; i++) {
	mesh->sol[k].m[i] = fbuf[off+i];
	mesh->sol[k].bb  += fbuf[off+i]*fbuf[off+i];
      }
      mesh->sol[k].bb = sqrt(mesh->sol[k].bb);
      if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
      if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
    }
    
    //printf("max= %f, min= %f",mesh->bbmin,mesh->bbmax);
    break;

  case GmfSymMat:
    if ( ddebug )  printf("   metric field\n");
    mesh->nfield = sol->dim*(sol->dim+1) / 2;
    
    
    for (k=1; k<=nel; k++) {
      // reading data must be double !!! 
      getline_bin_double_noref(sol->dim*(sol->dim+1)/2, TenSol);
     
      for (i=0; i<sol->dim*(sol->dim+1)/2; i++) {
	fbuf[off+i] = TenSol[i];
      }
      
      if ( sol->dim == 2 ) {
	for (i=0; i<3; i++)
	  mesh->sol[k].m[i] = m[i] = fbuf[off+i];
	iord = eigen2(m,lambda,vp);
	mesh->sol[k].bb = min(lambda[0],lambda[1]);
	if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
	if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
      }
      else {
	for (i=0; i<6; i++)
	  mesh->sol[k].m[i] = fbuf[off+i];
	mesh->sol[k].m[2] = fbuf[off+3];
	mesh->sol[k].m[3] = fbuf[off+2];
	for (i=0; i<6; i++)  m[i] = mesh->sol[k].m[i];
	iord = eigenv(1,m,lambda,eigv);
	if ( iord ) {
	  mesh->sol[k].bb = lambda[0];
	  mesh->sol[k].bb = max(mesh->sol[k].bb,lambda[1]);
	  mesh->sol[k].bb = max(mesh->sol[k].bb,lambda[2]);
	  if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
	  if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
	}
      }
    }
    
    break;
  }
  retcode=1;
  
  return retcode;
}

/*load solution (metric) */
int loadSol_popen_bin(pMesh mesh,char *filename,int numsol) {
  pSolution    sol;
  double       dbuf[ GmfMaxTyp ];
  float        fbuf[ GmfMaxTyp ];
  double       m[6],lambda[3],eigv[3][3],vp[2][2];
  int          inm,k,i,key,nel,size,type,iord,off,typtab[GmfMaxTyp],ver,dim;
  char        *ptr,data[128];

  // rajout pour popen
  int       NumberofSolAT;
  char       *natureread;
  // rajout binaire
  int       KwdCod;
  int       cod;
  int       NulPos;
  int retcode=0;
  NumberofSolAT=0;
#ifdef WIN32
  _setmode(fileno(stdin),O_BINARY);     
#endif
  // 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 *)&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 *)&dim ,WrdSiz, 1, stdin);
  natureread="Dimension";
  printf(".sol: %s %i (mesh)%i (lecture)%i \n",natureread,dim,mesh->dim,ver);
  /*control of the dimension*/
  if( dim != mesh->dim ){
    fprintf(stderr,"  %%%% Wrong dimension %d.\n",dim);
    retcode=0;
    goto Lret; 
    //    return(0);
  }

  while( !feof(stdin) ){
   
    fread( (unsigned char *)&KwdCod ,WrdSiz, 1, stdin);
   
    if(KwdCod == GmfSolAtVertices){   
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      fread( (unsigned char *)&nel ,WrdSiz, 1, stdin);
      natureread = "SolAtVertices";
      fprintf(stdout,"SolAtVertices : nel %i, mesh->np %i \n",nel,mesh->np);
      
      if ( nel != mesh->np ) {
	fprintf(stderr,"  %%%% Wrong number: %d Solutions discarded\n",nel-mesh->np);
	retcode=0;
	goto Lret; 
	//	return(0);
      }
      mesh->typage = 2;
      key = GmfSolAtVertices;

      /*  type,size,typetab  */  
      read_TypeSizeTyptab_bin( &type, &size, typtab);
      printf("sol: %s; type %i; size%i;\n",natureread, type, size); 
      fflush(stdout);
      /* Reading solutions*/
      loadScaVecTen_bin( mesh, 1, dim, ver, nel, type, size, typtab, key);
    }
    
    if( mesh->dim == 2 && mesh->nt ){
      if(KwdCod == GmfSolAtTriangles){
	natureread = "SolAtTriangles";
	fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
	fread( (unsigned char *)&nel ,WrdSiz, 1, stdin);
	fprintf(stdout,"SolAtTriangles : nel %i, mesh->nt %i \n",nel,mesh->nt);
	if ( nel && nel != mesh->nt ) {
	  fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
	  retcode=0;
	  goto Lret; 
	  // return(0);
        }

	mesh->typage = 1;
	key = GmfSolAtTriangles;

	/*  type,size,typetab  */  
	read_TypeSizeTyptab_bin( &type, &size, typtab);
	printf("sol: %s; type %i; size%i;\n",natureread, type, size); 

	/* Reading solutions*/
	loadScaVecTen_bin( mesh, 1, dim, ver, nel, type, size, typtab, key);

      }
      
    }

    if( mesh->dim == 2 && mesh->nq ){
      if( KwdCod == GmfSolAtQuadrilaterals ){	
	natureread = "SolAtQuadrilaterals";
	fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
	fread( (unsigned char *)&nel ,WrdSiz, 1, stdin);
	fprintf(stdout,"SolAtQuadrilaterals : nel %i, mesh->nq %i \n",nel,mesh->nq);
	if ( nel && nel != mesh->nq ) {
	  fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
	  retcode=0;
	  goto Lret; 
	  //return(0);
	}
	
	mesh->typage = 1;
	key = GmfSolAtQuadrilaterals;

	/*  type,size,typetab  */  
	read_TypeSizeTyptab_bin( &type, &size, typtab);

	/* Reading solutions*/
	loadScaVecTen_bin( mesh, 1, dim, ver, nel, type, size, typtab, key);
      }
    }

    if( mesh->dim == 3 && mesh->ntet ){
      if( KwdCod == GmfSolAtTetrahedra ){
	natureread = "SolAtTetrahedra";
	fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
	fread( (unsigned char *)&nel ,WrdSiz, 1, stdin);
	fprintf(stdout,"SolAtTetrahedra : nel %i, mesh->ntet %i \n",nel,mesh->ntet);
	if ( nel && nel != mesh->ntet ) {
	  fprintf(stderr,"  %%%% Wrong number %d.\n",nel); 
	  retcode=0;
	  goto Lret; 
	  //return(0);
	}
	mesh->typage = 1;
	key = GmfSolAtTetrahedra;

	/*  type,size,typetab  */  
	read_TypeSizeTyptab_bin( &type, &size, typtab);

	/* Reading solutions */
	loadScaVecTen_bin( mesh, 1, dim, ver, nel, type, size, typtab, key);

      }
    }

    if( mesh->dim == 3 && mesh->nhex ){
      if(  KwdCod == GmfSolAtHexahedra ){
	natureread = "SolAtHexahedra";
	fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
	fread( (unsigned char *)&nel ,WrdSiz, 1, stdin);
	fprintf(stdout,"SolAtHexahedra : nel %i, mesh->nhex %i  %s\n",nel,mesh->nhex);
	if ( nel && nel != mesh->nhex ) {
	  fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
	  GmfCloseMesh(inm);
	  retcode=0;
	  goto Lret; 
	  //return(0);
	}
	mesh->typage = 1;
	key = GmfSolAtHexahedra;

	/*  type,size,typetab  */  
	read_TypeSizeTyptab_bin( &type, &size, typtab);

	/* Reading solutions*/
	loadScaVecTen_bin( mesh, 1, dim, ver, nel, type, size, typtab, key);
      }
    }
    if( KwdCod == GmfEnd ){
      fread( (unsigned char *)&NulPos ,WrdSiz, 1, stdin);
      printf("End of solution\n");
      if( ddebug ) printf("Reading of mesh file is finished");
      break;
    }    
  }
  retcode=1;
 Lret:
#ifdef WIN32
  _setmode(fileno(stdin),O_BINARY);     
#endif 
  return(retcode);
}