medit.cpp

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

    vzy.init(datasize);
    vzz.init(datasize);
    
    if( nargs[2] ){
      boolsave= true;
      saveff = GetAny<string *>( (*nargs[2])(stack) );
      int ddim = 3;
    
      for (size_t i=0;i<offset;i++){
	solnbfloat = solnbfloat + l[i].nbfloat;
      }
      
      solsave.init(solnbfloat,datasize);
      //cout << "solsave.size()= " << solsave.size()  << endl;
      solsave=0.;
      //cout << solsave  << endl;
    }

  }

  int nboftmp = 0;
    
  FILE *popenstream= popen(commandline,MODE_WRITE_BINARY);
  if( !popenstream){
    cerr << " Error popen : " << commandline<<endl;
    exit(1);
  }
 
  // mesh
  int jojo1;
  for(int jojo=0; jojo<smedit; jojo++){

    if(filebin){
      int cod = 1;
      int KwdCod;
      int NulPos = 0;
      
      // determination of number solutions associated with a mesh
      fwrite( (unsigned char *) &cod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &ver, WrdSiz,1,popenstream );
      KwdCod = GmfDimension;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &dimp, WrdSiz,1,popenstream );
     

      // vertex 
      KwdCod = GmfVertices;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &nv, WrdSiz,1,popenstream );   
      for (int k=0; k<nv; k++) {
	const  Vertex3 & P = Th.vertices[k];
	fx=P.x; fy=P.y; fz=P.z;
	fwrite( (unsigned char *) &fx, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &fy, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &fz, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &(P.lab), WrdSiz,1,popenstream );
	//fprintf(popenstream,"%f %f %i\n",fx,fy,P.lab);
      }
    
      // tetrahedra
      KwdCod= GmfTetrahedra;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &nt, WrdSiz,1,popenstream );  
      for (int k=0; k<nt; k++) {
	const Tet & K(Th.elements[k]);
	int i0=Th.operator()(K[0])+1;
	int i1=Th.operator()(K[1])+1;
	int i2=Th.operator()(K[2])+1;
	int i3=Th.operator()(K[3])+1;
	int lab=K.lab;
	fwrite( (unsigned char *) &i0, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &i1, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &i2, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &i3, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &lab, WrdSiz,1,popenstream );
	//fprintf(popenstream,"%i %i %i %i %i\n",i0,i1,i2,i3,lab);
      }

      // triangles
      KwdCod = GmfTriangles;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &nbe, WrdSiz,1,popenstream );
      for (int k=0; k<nbe; k++) {
	const Triangle3 & K(Th.be(k));
	int i0=Th.operator()(K[0])+1;
	int i1=Th.operator()(K[1])+1;
	int i2=Th.operator()(K[2])+1;
	int lab=K.lab;
	fwrite( (unsigned char *) &i0, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &i1, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &i2, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &lab, WrdSiz,1,popenstream );
	//fprintf(popenstream,"%i %i %i %i\n",i0,i1,i2,lab);
      }
            
      // End
      KwdCod = GmfEnd;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      //fprintf(popenstream,"End\n");
    }
    else{
      fprintf(popenstream,"MeshVersionFormatted\n");
      fprintf(popenstream,"%i\n",ver);
      fprintf(popenstream,"Dimension\n");
      fprintf(popenstream,"%i\n",dimp);
      fprintf(popenstream,"Vertices\n");
      fprintf(popenstream,"%i\n",nv);
      
      for (int k=0; k<nv; k++) {
	const  Vertex3 & P = Th.vertices[k];
	fx=P.x; fy=P.y; fz=P.z;
	fprintf(popenstream,"%f %f %f %i\n",fx,fy,fz,P.lab);
      }
      fprintf(popenstream,"Tetrahedra\n");
      fprintf(popenstream,"%i\n",nt);
      for (int k=0; k<nt; k++) {
	const Tet & K(Th.elements[k]);
	int i0=Th.operator()(K[0])+1;
	int i1=Th.operator()(K[1])+1;
	int i2=Th.operator()(K[2])+1;
	int i3=Th.operator()(K[3])+1;
	int lab=K.lab;
	fprintf(popenstream,"%i %i %i %i %i\n",i0,i1,i2,i3,lab);
      }
      fprintf(popenstream,"Triangles\n");
      fprintf(popenstream,"%i\n",nbe);
      for (int k=0; k<nbe; k++) {
	const Triangle3 & K(Th.be(k));
	int i0=Th.operator()(K[0])+1;
	int i1=Th.operator()(K[1])+1;
	int i2=Th.operator()(K[2])+1;
	int lab=K.lab;
	fprintf(popenstream,"%i %i %i %i\n",i0,i1,i2,lab);
      }
      fprintf(popenstream,"End");
    }
    
    if(nbsol > 0){
      if(filebin){
	int cod = 1;
	int NulPos = 0;
	int KwdCod;
	fwrite( (unsigned char *) &cod, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &ver, WrdSiz,1,popenstream );
	KwdCod = GmfDimension;
	fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &dimp, WrdSiz,1,popenstream );
      }
      else{
	fprintf(popenstream,"MeshVersionFormatted %i\n",ver);
	fprintf(popenstream,"Dimension %i\n",dimp);
	
	// detemination of solution
	// default scalaire // faire tableau pour plusieurs
      }
      
      typsol = TypTab[jojo];
	
      //fprintf(popenstream,"%i %i\n",1,typsol);
      
      
      if(order == 0){
	if(filebin){
	  int NulPos = 0;
	  int KwdCod = GmfSolAtTetrahedra;
	  int codtypjm = 1;
	  fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &datasize, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &codtypjm, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &typsol, WrdSiz,1,popenstream );
	}
	else{
	  fprintf(popenstream,"SolAtTetrahedra\n");
	  fprintf(popenstream,"%i\n",datasize);
	  fprintf(popenstream,"%i %i\n",1,typsol);
	}
	if(typsol==1){
	  //KN<double> solsca(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  R3 Cdg_hat = R3(1./4.,1./4.,1./4.);  

	  vxx=0.;
	  for (int it=0;it<nt;it++){
	    jojo1 = jojo+1+offset*numTht[it];
	    const Tet & K(Th.elements[it]);
	    mp3->set( Th, K(Cdg_hat), Cdg_hat, K, K.lab);
	    
	    vxx[it] = l[jojo1].eval(0,stack); //GetAny< double >( (*nargs[1])(stack) );
	    
	  }
	  
	  if(filebin){
	    for(int k=0; k<nt; k++){
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz,2,popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nt; k++){
	      fprintf(popenstream,"%f\n",vxx[k]);
	    }
	  }
	}
	else if(typsol==2){
	  //KN<double> vxx(nv),vyy(nv),vzz(nv);
	  vxx=0.; vyy=0.; vzz=0.;    
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  R3 Cdg_hat = R3(1./4.,1./4.,1./4.); 
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    const Tet & K(Th.elements[it]);
	    mp3->set( Th, K(Cdg_hat), Cdg_hat, K, K.lab);
	    
	    vxx[it] = l[jojo1].eval(0,stack); //GetAny< double >( (*xx)(stack) );
	    vyy[it] = l[jojo1].eval(1,stack); //GetAny< double >( (*yy)(stack) );
	    vzz[it] = l[jojo1].eval(2,stack); //GetAny< double >( (*zz)(stack) );
	    
	  }
	  if(filebin){
	    for(int k=0; k<nt; k++){
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vyy[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vzz[k]), WrdSiz,2,popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nt; k++){
	      fprintf(popenstream,"%f %f %f\n", vxx[k], vyy[k], vzz[k]);
	    }
	  }
	}
	else if(typsol==3){
	  //KN<double> vxx(nv),vyx(nv),vyy(nv),vzx(nv),vzy(nv),vzz(nv);
	  vxx(nv)=0.;vyx(nv)=0.;vyy(nv)=0.;
	  vzx(nv)=0.;vzy(nv)=0.;vzz(nv)=0.;
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  R3 Cdg_hat = R3(1./4.,1./4.,1./4.); 
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    const Tet & K(Th.elements[it]);
	    mp3->set( Th, K(Cdg_hat), Cdg_hat, K, K.lab);
	    
	    vxx[it] = l[jojo1].eval(0,stack); //GetAny< double >( (*tsxx)(stack) );
	    vyx[it] = l[jojo1].eval(1,stack); //GetAny< double >( (*tsyx)(stack) );
	    vyy[it] = l[jojo1].eval(2,stack); //GetAny< double >( (*tsyy)(stack) );
	    vzx[it] = l[jojo1].eval(3,stack); //GetAny< double >( (*tszx)(stack) );
	    vzy[it] = l[jojo1].eval(4,stack); //GetAny< double >( (*tszy)(stack) );
	    vzz[it] = l[jojo1].eval(5,stack); //GetAny< double >( (*tszz)(stack) );
	    
	  }
	  if(filebin){
	    for(int k=0; k<nt; k++){
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vyx[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vyy[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vzx[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vzy[k]), WrdSiz,2,popenstream );
	      fwrite( (unsigned char *) &(vzz[k]), WrdSiz,2,popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nt; k++){
	      fprintf(popenstream,"%f %f %f %f %f %f\n",vxx[k],vyx[k],vyy[k],vzx[k],vzy[k],vzz[k]);
	    } 
	  }
	}
      }
      else if(order == 1){
	if(filebin){
	  int NulPos = 0;
	  int KwdCod = GmfSolAtVertices;
	  int codtypjm = 1;
	  fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &datasize, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &codtypjm, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &typsol, WrdSiz,1,popenstream );
	}
	else{
	  fprintf(popenstream,"SolAtVertices\n");
	  fprintf(popenstream,"%i\n",datasize);
	  fprintf(popenstream,"%i %i\n",1,typsol);
	}
	if(typsol==1){
	  //KN<double> solsca(nv);
	  KN<int> takemesh(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  	  
	  takemesh=0;
	  vxx=0.;

	  for (int it=0;it<Th.nt;it++){
	    jojo1 = jojo+1+offset*numTht[it];
	    for( int iv=0; iv<4;iv++){
	      int i=Th(it,iv);  
	      
	      mp3->setP(&Th,it,iv);
	      vxx[i] = vxx[i]+l[jojo1].eval(0,stack); //GetAny< double >( (*nargs[1])(stack) );
	      
	      takemesh[i] = takemesh[i]+1;
	    }
	  }
	  if(filebin){
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];	      
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz,2,popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];
	      fprintf(popenstream,"%f\n",vxx[k]);
	    }
	  }
	}
	else if(typsol==2){
	  //KN<double> vxx(nv),vyy(nv),vzz(nv);
	  KN<int> takemesh(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  
	  vxx=0.;
	  vyy=0.;
	  vzz=0.;
	  takemesh=0;
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    for( int iv=0; iv<4; ++iv){
	      int i=Th(it,iv);  
	      
	      //if(takemesh[i]==0){
	      mp3->setP(&Th,it,iv);
	      vxx[i] = vxx[i]+l[jojo1].eval(0,stack); //GetAny< double >( (*xx)(stack) );
	      vyy[i] = vyy[i]+l[jojo1].eval(1,stack); //GetAny< double >( (*yy)(stack) );
	      vzz[i] = vzz[i]+l[jojo1].eval(2,stack); //GetAny< double >( (*zz)(stack) );
	      
	      takemesh[i] = takemesh[i]+1;
	      //}
	    }
	  }
	  if(filebin){
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];
	      vyy[k]=vyy[k]/takemesh[k];
	      vzz[k]=vzz[k]/takemesh[k];
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vyy[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vzz[k]), WrdSiz, 2, popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];
	      vyy[k]=vyy[k]/takemesh[k];
	      vzz[k]=vzz[k]/takemesh[k];
	      fprintf(popenstream,"%f %f %f\n", vxx[k], vyy[k], vzz[k]);
	    }
	  }
	}
	else if(typsol==3){
	  //KN<double> vxx(nv),vyx(nv),vyy(nv),vzx(nv),vzy(nv),vzz(nv);
	  KN<int> takemesh(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 

	  vxx=0.;
	  vyx=0.;
	  vyy=0.;
	  vzx=0.;
	  vzy=0.;
	  vzz=0.;
	  takemesh=0;
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    for( int iv=0; iv<4; ++iv){
	      int i=Th(it,iv);  
	      
	      mp3->setP(&Th,it,iv);
	      vxx[i] = vxx[i] + l[jojo1].eval(0,stack); //GetAny< double >( (*tsxx)(stack) );
	      vyx[i] = vyx[i] + l[jojo1].eval(1,stack); //GetAny< double >( (*tsyx)(stack) );
	      vyy[i] = vyy[i] + l[jojo1].eval(2,stack); //GetAny< double >( (*tsyy)(stack) );
	      vzx[i] = vzx[i] + l[jojo1].eval(3,stack); //GetAny< double >( (*tszx)(stack) );
	      vzy[i] = vzy[i] + l[jojo1].eval(4,stack); //GetAny< double >( (*tszy)(stack) );
	      vzz[i] = vzz[i] + l[jojo1].eval(5,stack); //GetAny< double >( (*tszz)(stack) );
	      
	      takemesh[i] = takemesh[i]+1; 	
	    }
	  }
	  if(filebin){
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];
	      vyx[k]=vyx[k]/takemesh[k];
	      vyy[k]=vyy[k]/takemesh[k];
	      vzx[k]=vzx[k]/takemesh[k];
	      vzy[k]=vzy[k]/takemesh[k];
	      vzz[k]=vzz[k]/takemesh[k];

	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vyx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vyy[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vzx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vzy[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vzz[k]), WrdSiz, 2, popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];
	      vyx[k]=vyx[k]/takemesh[k];
	      vyy[k]=vyy[k]/takemesh[k];
	      vzx[k]=vzx[k]/takemesh[k];
	      vzy[k]=vzy[k]/takemesh[k];
	      vzz[k]=vzz[k]/takemesh[k];
	      
	      fprintf(popenstream,"%f %f %f %f %f %f\n",vxx[k], vyx[k], vyy[k], vzx[k], vzy[k], vzz[k]);
	    }
	  }

	}
      }

      if(filebin){
	int NulPos = 0;
	int KwdCod = GmfEnd;
	fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      }
      else{
	fprintf(popenstream,"End");
      }



      if(boolsave){
	if(verbosity) cout << "writing solution in file" << endl;
	if(typsol==1){
	  writetabsol( datasize, nboftmp, vxx, solsave);
	  nboftmp=nboftmp+1;
	}
	else if(typsol==2){
	  writetabsol( datasize, nboftmp, vxx, vyy, vzz, solsave);
	  nboftmp=nboftmp+3;
	}
	else if(typsol==3){
	  writetabsol( datasize, nboftmp, vxx, vyx, vyy, vzx, vzy, vzz, solsave);
	  nboftmp=nboftmp+6;
	}
	//cout << "finish writing solution in file" << endl;
	//cout << "datasize=" << datasize << " " << "solnbfloat=" << solnbfloat << endl;
      }

    }
  }
  // fermeture du stream pour popen
  bool wait=TheWait;
  if (nargs[3]) wait= GetAny<bool>((*nargs[3])(stack));
  
  bool plotting = true;   
  //  drawing part  ------------------------------
  
  if (wait && !NoWait) 
    {
      pclose(popenstream);
    }
  else
    {
      fclose(popenstream);
    }
  
  if(boolsave){
    int outm;
    int nbtype=nbsol;
    float *OutSolTab = new float[solnbfloat];
    
    if ( !(outm = GmfOpenMesh(saveff->c_str(),GmfWrite,ver,3)) ) {
      cerr <<"  -- Mesh3::Save  UNABLE TO OPEN  :"<< saveff << endl;
      exit(1);
    }
    
    if(order == 0){	  
      GmfSetKwd(outm,GmfSolAtTetrahedra, datasize, nbtype, TypTab);
      for (int k=0; k<datasize; k++){
	for (int i=0; i<solnbfloat ;i++){
	  OutSolTab[i] =  solsave(i,k);
	}
	GmfSetLin(outm, GmfSolAtTetrahedra, OutSolTab); 
      }
    }
    else if(order == 1){
      GmfSetKwd(outm,GmfSolAtVertices, datasize, nbtype, TypTab);
      for (int k=0; k<datasize; k++){
	for (int i=0; i<solnbfloat ;i++){
	  OutSolTab[i] =  solsave(i,k);
	}
	GmfSetLin(outm, GmfSolAtVertices, OutSolTab);
      }
    }
    delete [] OutSolTab;
  }
  return valsortie;
}

//  truc pour que la fonction 
// Init::Init() soit appele a moment du chargement dynamique
// du fichier 
//  

class Init { public:
  Init();
};

static Init init;  //  une variable globale qui serat construite  au chargement dynamique 

Init::Init(){  // le constructeur qui ajoute la fonction "splitmesh3"  a freefem++ 
  
  typedef Mesh *pmesh;
  //typedef Mesh2 *pmesh2;
  typedef Mesh3 *pmesh3;
  
  if (verbosity)
    cout << " load:popen.cpp  " << endl;
  
  // 2D
  Global.Add("medit","(",new OneOperatorCode<PopenMeditMesh_Op>);
  Global.Add("savesol","(",new OneOperatorCode<datasolMesh2_Op> );
  
  // 3D
  Global.Add("medit","(",new OneOperatorCode< PopenMeditMesh3_Op<v_fes3> >);
  Global.Add("savesol","(",new OneOperatorCode< datasolMesh3_Op<v_fes3> >);
  
}