medit.cpp

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

	typsol = TypTab[jojo];
	
	fprintf(popenstream,"%i %i\n",1,typsol);
      }

      
      if(typsol==1){
	if(order==0){
	  
	  vxx=0.;
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  R2 Cdg_hat = R2(1./3.,1./3.);  
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    const Mesh::Triangle & K(Th.t(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(order==1){
	  //KN<double> solsca(nv);    
	  vxx=0.;
	  KN<int> takemesh(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  
	  takemesh=0;
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    for( int iv=0; iv<3; ++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){
	if(order==0){	
	  //KN<double>  vxx(nt),vyy(nt);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  R2 Cdg_hat = R2(1./3.,1./3.);  
	  
	  vxx=0.; vyy=0.;
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    const Mesh::Triangle & K(Th.t(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) );
	  }
	  
	  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 );
	    }
	  }
	  else{
	    for(int k=0; k<nt; k++){
	      fprintf(popenstream,"%f %f\n",vxx[k],vyy[k]);
	    }
	  }
	}
	
	else if(order==1){
	  //KN<double> vxx(nv),vyy(nv);
	  KN<int> takemesh(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  
	  takemesh=0;
	  vxx=0.; vyy=0.;
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    for( int iv=0; iv<3; ++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) );
	      
	      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];
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vyy[k]), WrdSiz, 2, popenstream );
	    }
	  }
	  else{
	    for(int k=0; k<nv; k++){
	      vxx[k]=vxx[k]/takemesh[k];
	      vyy[k]=vyy[k]/takemesh[k];
	      fprintf(popenstream,"%f %f\n", vxx[k], vyy[k]);
	    }
	  }
	  
	}
      }
      else if(typsol==3){
	if(order==0){
	  //KN<double>  vxx(nt),vyx(nt),vyy(nt);
	  vxx=0.; vyx=0.; vyy=0.;
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  R2 Cdg_hat = R2(1./3.,1./3.);  
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    const Mesh::Triangle & K(Th.t(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) );
	  }
	  
	  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 );
	    }
	  }
	  else{
	    for(int k=0; k<nt; k++){
	      fprintf(popenstream,"%f %f %f\n", vxx[k], vyx[k], vyy[k]);
	    }
	  }
	}
	else if(order==1){
	  //KN<double> vxx(nv),vyx(nv),vyy(nv);
	  KN<int> takemesh(nv);
	  MeshPoint *mp3(MeshPointStack(stack)); 
	  
	  takemesh=0;
	  vxx=0.; vyx=0.; vyy=0.;
	  
	  for (int it=0;it<Th.nt;++it){
	    jojo1 = jojo+1+offset*numTht[it];
	    for( int iv=0; iv<3; ++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 >( (*tsxx)(stack) );
	      vyx[i] = vyx[i]+l[jojo1].eval(0,stack); //GetAny< double >( (*tsyx)(stack) );
	      vyy[i] = vyy[i]+l[jojo1].eval(0,stack);  //GetAny< double >( (*tsyy)(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];
	      
	      fwrite( (unsigned char *) &(vxx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vyx[k]), WrdSiz, 2, popenstream );
	      fwrite( (unsigned char *) &(vyy[k]), WrdSiz, 2, popenstream );
	      
	      //fprintf(popenstream,"%f %f %f\n",vxx[k],vyx[k],vyy[k]);
	    }
	  }
	  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];
	      fprintf(popenstream,"%f %f %f\n",vxx[k],vyx[k],vyy[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\n");
      }

      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, solsave);
	  nboftmp=nboftmp+2;
	}
	else if(typsol==3){
	  writetabsol( datasize, nboftmp, vxx, vyx, vyy, solsave);
	  nboftmp=nboftmp+3;
	}
	/*cout << "fin writing solution in file nboftmp=" << nboftmp << endl;
	cout << "datasize=" << datasize << " " << "solnbfloat=" << solnbfloat << " boolsave=" << boolsave << endl;
	for(int iy=0; iy<datasize; iy++){
	  if(iy==0 || iy==datasize-1) cout << iy <<" " <<solsave(0,iy) << endl;	      
	}
	*/
      }
    }
  }

  // rajout pour wait
  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);
    }
  
  // rajout pout la sauvegarde de la solution
  if(boolsave){
    //cout <<"save solution in file avec printf\n"<< endl;;
    int outm;
    int nbtype=nbsol;
    float *OutSolTab = new float[solnbfloat];
    
    if ( !(outm = GmfOpenMesh(saveff->c_str(),GmfWrite,ver,2)) ) {
      cerr <<"  -- Mesh3::Save  UNABLE TO OPEN  :"<< saveff << endl;
      exit(1);
    }
    
    if(order == 0){	  
      GmfSetKwd(outm,GmfSolAtTriangles, datasize, nbtype, TypTab);
      for (int k=0; k<datasize; k++){
	for (int i=0; i<solnbfloat ;i++){
	  OutSolTab[i] =  solsave(i,k);
	}
	GmfSetLin(outm, GmfSolAtTriangles, 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;
  }
  
  delete pTh;

 
  return valsortie;
}

template<class v_fes>
class PopenMeditMesh3_Op : public E_F0mps 
{
public:
  typedef long  Result;
  Expression eTh;
  Expression filename;
  long offset;
  long nbTh;
  struct Expression2 {
    long what; // 0 mesh, 1 scalar, 2 vector, 3 symtensor
    long nbfloat; // 0 mesh(3D), 1 scalar, 2 vector (3D), 3 symtensor(3D)
    Expression e[6];
    Expression2() {e[0]=0; e[1]=0; e[2]=0;  e[3]=0; e[4]=0; e[5]=0; what=0; nbfloat=0;};
    Expression &operator[](int i){return e[i];}
    double eval(int i,Stack stack) const  { 
    if (e[i]) {
      return GetAny< double >( (*e[i])(stack) );
    }
    else 
      return 0;
    }
    const Mesh3 & evalm(int i,Stack stack) const  { throwassert(e[i]);return  * GetAny< pmesh3 >((*e[i])(stack)) ;}
  };
  vector<Expression2> l;

  static const int n_name_param = 5;  
  static basicAC_F0::name_and_type name_param[] ;
  Expression nargs[n_name_param];
  long arg(int i,Stack stack,int a) const{ return nargs[i] ? GetAny< long >( (*nargs[i])(stack) ): a;}
  string*  arg(int i,Stack stack, string* a ) const{ return nargs[i] ? GetAny< string* >( (*nargs[i])(stack) ): a;}

public:
  PopenMeditMesh3_Op(const basicAC_F0 &  args) : l( args.size()-1 )
  {
    int nbofsol;
    int ddim=3;
    int stsize=6;
    
    args.SetNameParam(n_name_param,name_param,nargs);   
   
    if (BCastTo<string *>(args[0])) filename = CastTo<string *>(args[0]);
    //if (BCastTo<pmesh3>(args[1])) eTh= CastTo<pmesh3>(args[1]);
    
    for (size_t i=1;i<args.size();i++){
      size_t jj=i-1;
      if ( BCastTo<double>(args[i]) )
	{
	  l[jj].what=1;
	  l[jj].nbfloat=1;
	  l[jj][0]=to<double>( args[i] );
	}
      else if ( args[i].left()==atype<E_Array>() )
	{
	  const E_Array * a0  = dynamic_cast<const E_Array *>( args[i].LeftValue() );
	  //cout << "taille" << a0->size() << endl;
	  //if (a0->size() != ddim || a0->size() != stsize) 
	  //  CompileError("savesol in 2D: vector solution is 2 composant, vector solution is 3 composant");

	  if( a0->size() == ddim){
	    // vector solution
	    l[jj].what=2;
	    l[jj].nbfloat=ddim;
	    for(int j=0; j<ddim; j++){
	      l[jj][j] = to<double>( (*a0)[j]);
	    }
	  }
	  else if( a0->size() == stsize){
	    // symmetric tensor solution
	    l[jj].what=3;
	    l[jj].nbfloat=stsize;
	    for(int j=0; j<stsize; j++){
	      l[jj][j] = to<double>( (*a0)[j]);
	    }
	  }
	  
	}
      else if( BCastTo<pmesh3>(args[i]) ){
	l[jj].what    = 0;
	l[jj].nbfloat = 0;
	l[jj][0] = CastTo<pmesh3>(args[i]);
      }
      else {
	CompileError("medit 3d: Sorry no way to save this kind of data");
      }      
    }
    // determination of the number of solutions.
    size_t lastTh=0;
    long offset1;
    offset=0;
    nbTh=0;
    for(size_t jj=0; jj<l.size(); jj++){
      if(l[jj].what==0){
	nbTh++;
	offset1=jj-lastTh;
	if(offset==0){
	  offset=offset1;
	}
	else if(offset != offset1){
	  CompileError("the number of solution by mesh is different");
	}
      }
    }
    if(offset==0) offset=l.size(); 
    
  }
  
  static ArrayOfaType  typeargs() { return  ArrayOfaType( atype<string *>(), atype<pmesh3>(), true); }// all type
  static  E_F0 * f(const basicAC_F0 & args) { return new PopenMeditMesh3_Op(args);} 
  AnyType operator()(Stack stack)  const ;
};
template<class v_fes>
basicAC_F0::name_and_type PopenMeditMesh3_Op<v_fes>::name_param[]= {
  {  "order", &typeid(long)},
  {  "meditff", &typeid(string*)},
  {  "save",&typeid(string*)},
  {  "wait",&typeid(bool)},
  {  "bin",&typeid(long)}
};

template<class v_fes>
AnyType PopenMeditMesh3_Op<v_fes>::operator()(Stack stack)  const 
{
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  long order (arg(0,stack,1));
  //
  int ver = GmfFloat;
  int dimp =3;
  float fx,fy,fz;
  //
  long valsortie=0;
  int typsol,nbsol;
  nbsol = offset-1;
  
  int TypTab[l.size()-1];
  for (size_t i=0;i<l.size()-1;i++){
    TypTab[i]=l[i+1].what;
  }

  //  string stringemptymedit= string("medit");
  string * ffname  = GetAny<string *>( (*filename)(stack) );
  string * meditff(arg(1,stack,&stringffmedit));

  long filebin (arg(4,stack,1));
  int smedit=max(1,nbsol);     
  char * commandline = meditcmd( filebin, nbsol, smedit, *meditff, *ffname);
  
  printf("version de medit %s\n",commandline);  
  if(verbosity) cout << "number of solution = " << offset-1 << endl;
  if(verbosity) cout << "number of mesh     = " << nbTh << endl;

  // lecture des differents maillages
  int nv=0,nt=0,nbe=0; // sommet, triangles, arretes du maillage unifies

  //cout << "commencement du maillage " << endl;
  for(size_t i=0; i<l.size();i=i+offset){
    if(l[i].what!=0)
      cerr << "this element is not a mesh" << i << endl;
    const Mesh3 &Thtmp =l[i].evalm(0,stack);
    nv  += Thtmp.nv;
    nt  += Thtmp.nt;
    nbe += Thtmp.nbe;  
    //cout << "valeur de i=" << i << "l.size()=" << l.size() << endl;
    //cout << "vertex "<< nv << " tetrahedra "<< nt << " triangle " << nbe << endl;  
  }
  
  Vertex3    *v = new Vertex3[nv];
  Tet        *t = new Tet[nt];
  Triangle3  *b = new Triangle3[nbe]; 
  Tet       *tt = t;
  Triangle3 *bb = b;

  int iv=0,it=0,ibe=0;
  int numTht[nt]; // numero of Th assoctiated with a tetrahedra 

  int jt=0;
  for(size_t i=0; i<l.size();i=i+offset){
    int nvtmp  = iv;
    int nttmp  = it;
    int nbetmp = ibe;
      
    const Mesh3 &Thtmp =l[i].evalm(0,stack);
    for (int ii=0;ii<Thtmp.nv;ii++){
      const Vertex3 &vi(Thtmp.vertices[ii]);
      v[iv].x = vi.x;
      v[iv].y = vi.y;
      v[iv].z = vi.z;
      v[iv].lab = vi.lab;
      iv++;
    }

    for (int ii=0;ii<Thtmp.nt;ii++){
      const Tet &vi(Thtmp.elements[ii]);
      int iv[4];
      iv[0] = nvtmp + Thtmp.operator()(vi[0]);
      iv[1] = nvtmp + Thtmp.operator()(vi[1]);
      iv[2] = nvtmp + Thtmp.operator()(vi[2]);
      iv[3] = nvtmp + Thtmp.operator()(vi[3]);
      (*tt++).set(v,iv,vi.lab);
      numTht[it] = jt;
      it++;      
    }

    for (int ii=0;ii<Thtmp.nbe;ii++){
      const  Triangle3 &vi(Thtmp.be(ii));
      int iv[3];
      iv[0] = nvtmp + Thtmp.operator()(vi[0]);
      iv[1] = nvtmp + Thtmp.operator()(vi[1]);
      iv[2] = nvtmp + Thtmp.operator()(vi[2]);
      (*bb++).set(v,iv,vi.lab);
      ibe++;
    }
    jt++;
  }
  assert( it==nt ); assert(iv==nv); assert(ibe=nbe);
  if(verbosity) cout << "meditff :: Value of elements: vertex "<< nv << " Tet "<< nt << " triangle " << nbe << endl;  

  Mesh3 *pTh = new Mesh3(nv,nt,nbe,v,t,b);
  Mesh3 &Th = *pTh;

  //cout << "Mesh is created" << endl;
  // determination of the number of elements to represent the solution
  int datasize;
  if(order == 0) datasize= nt;
  if(order == 1) datasize= nv;

  // cas de sauvegarde
  bool boolsave = false;
  int solnbfloat=0;
  KNM<double> solsave(1,1);
  string * saveff;
  KN<double> vxx,vyx,vyy,vzx,vzy,vzz;
 
  //cout << "nbsol=" << endl;

  if(nbsol > 0){
    
    vxx.init(datasize);
    vyx.init(datasize);
    vyy.init(datasize);
    vzx.init(datasize);