medit.cpp

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

	    for(size_t ii=0;ii<l.size();ii++){
	      for(size_t j=0;j<l[ii].nbfloat;j++){
		//cout << "ii=" << ii << " j=" << j<< endl;
		valsol[i*solnbfloat+h] = l[ii].eval(j,stack);
		h=h+1;
	      }
	    } 
	    assert(solnbfloat==h);
	    takemesh[i] = takemesh[i]+1;
	}
      }
    }


    GmfSetKwd(outm,GmfSolAtVertices, nbsol, nbtype, TypTab);
    for (int k=0; k<nbsol; k++){
      for (int i=0; i<solnbfloat ;i++){
	OutSolTab[i] =  valsol(k*solnbfloat+i);
      }
      GmfSetLin(outm, GmfSolAtVertices, OutSolTab);
    }
  }
  GmfCloseMesh(outm);
  delete [] OutSolTab;
  return longdefault;
}

//*************************
//
//  medit
//
//*************************

static char * meditcmd(long filebin, int nbsol, int smedit, const string &meditff, const string & ffnn)
{
  string meditcmm=meditff;
  int ddebug=0;
  if(ddebug)
    {
      meditcmm += ' ';
      meditcmm += medit_debug;
    }
  meditcmm += ' ';
  meditcmm += medit_popen;
  if(filebin)
    {
      meditcmm += ' ';
      meditcmm += medit_bin;
    }
  if(nbsol) 
    {
      meditcmm += ' ';
      meditcmm += medit_addsol;
    }
  
  char meditsol[5];
  sprintf(meditsol," %i",smedit);
  meditcmm += meditsol; 
   
  meditcmm += ' ';
  meditcmm += ffnn;
	  
  char * ret= new char[meditcmm.size()+1];
  //char ret[meditcmm.size()+1];
  strcpy( ret, meditcmm.c_str()); 
  return ret;
}

void writetabsol(const int &tsize, const int &nbofsol,const KN<double> &v1, KNM<double> &vv){
  
  for(int i=0; i<tsize; i++){
    vv(nbofsol,i)   = v1(i);
  }
}

void writetabsol(const int &tsize, const int &nbofsol,const KN<double> &v1, const KN<double> &v2, KNM<double> &vv){

  for(int i=0; i<tsize; i++){
    vv(nbofsol,i)   = v1(i);
    vv(nbofsol+1,i) = v2(i);
  }
}

void writetabsol(const int &tsize, const int &nbofsol,const KN<double> &v1, const KN<double> &v2, const KN<double> &v3, KNM<double> &vv){
  
  for(int i=0; i<tsize; i++){
    vv(nbofsol,i)   = v1(i);
    vv(nbofsol+1,i) = v2(i);
    vv(nbofsol+2,i) = v3(i);
  }

}

void writetabsol(const int &tsize, const int &nbofsol,const KN<double> &v1, const KN<double> &v2, const KN<double> &v3, 
	     const KN<double> &v4, const KN<double> &v5, const KN<double> &v6,KNM<double> &vv){
  
  for(int i=0; i<tsize; i++){
    vv(nbofsol,i)   = v1(i);
    vv(nbofsol+1,i) = v2(i);
    vv(nbofsol+2,i) = v3(i);
    vv(nbofsol+3,i) = v4(i);
    vv(nbofsol+4,i) = v5(i);
    vv(nbofsol+5,i) = v6(i);
  }

}

class PopenMeditMesh_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; // 1 scalar, 2 vector (2D), 3 symtensor(2D) 
    Expression e[3];
    Expression2() {e[0]=0; e[1]=0; e[2]=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 Mesh & evalm(int i,Stack stack) const  { throwassert(e[i]);return  * GetAny< pmesh >((*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:
  PopenMeditMesh_Op(const basicAC_F0 &  args) : l( args.size()-1 )
  {
    int nbofsol;
    int ddim=2;
    int stsize=3;
    
    args.SetNameParam(n_name_param,name_param,nargs);   
   
    if (BCastTo<string *>(args[0])) filename = CastTo<string *>(args[0]);
  
    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() );
	 
	  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<pmesh>(args[i]) ){
	l[jj].what    = 0;
	l[jj].nbfloat = 0;
	l[jj][0] = CastTo<pmesh>(args[i]);
      }
      else {
	CompileError("savesol in 2D: 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(); 
    // The number of solution is exactly:  offset-1

    // verification que la nature des solutions sont identiques pour les differents maillages.

    //cout << "number of solution = " << offset-1 << endl;
    //cout << "number of mesh     = " << nbTh << endl;
  } 
  
  static ArrayOfaType  typeargs() { return  ArrayOfaType( atype<string *>(), atype<pmesh>(), true); }// all type
  static  E_F0 * f(const basicAC_F0 & args) { return new PopenMeditMesh_Op(args);} 
  AnyType operator()(Stack stack)  const ;
};

basicAC_F0::name_and_type PopenMeditMesh_Op::name_param[]= {
  {  "order", &typeid(long)},
  {  "meditff", &typeid(string*)},
  {  "save",&typeid(string*)},
  {  "wait",&typeid(bool)},
  {  "bin",&typeid(long)}
};


AnyType PopenMeditMesh_Op::operator()(Stack stack)  const 
{
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  long order (arg(0,stack,1));
  //
  int ver = GmfFloat;
  int dimp =2;
  float fx,fy;
  //
  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 stringffmedit= string("medit.exe"); 
  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);

 
  // 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 Mesh &Thtmp =l[i].evalm(0,stack);
    nv  += Thtmp.nv;
    nt  += Thtmp.nt;
    nbe += Thtmp.neb;  
    //cout << "valeur de i=" << i << "l.size()=" << l.size() << endl;
    //cout << "vertex "<< nv << " triangle "<< nt << " edge " << nbe << endl;  
  }
  
  Mesh::Vertex        *v= new Mesh::Vertex[nv];
  Mesh::Triangle      *t= new Mesh::Triangle[nt];
  Mesh::BorderElement *b= new Mesh::BorderElement[nbe]; 
  Mesh::Triangle      *tt=t;
  Mesh::BorderElement *bb=b;
 
  int iv,it,ibe;
  iv=0;
  it=0;
  ibe=0;

  int numTht[nt]; // numero of Th assoctiated with a triangles 
  //int numTht[nbe]; // numero of Th assoctiated with a BoundaryEdge

  int jt=0;


  for(size_t i=0; i<l.size();i=i+offset){
    int nvtmp  = iv;
    int nttmp  = it;
    int nbetmp = ibe;
      
    const Mesh &Thtmp =l[i].evalm(0,stack);
    for (int ii=0;ii<Thtmp.nv;ii++){
      const Mesh::Vertex &vi(Thtmp(ii));
      v[iv] = vi;
      iv++;
    }

    for (int ii=0;ii<Thtmp.nt;ii++){
      const Mesh::Triangle &vi(Thtmp.t(ii));
      int i0 = nvtmp + Thtmp(ii,0);
      int i1 = nvtmp + Thtmp(ii,1);
      int i2 = nvtmp + Thtmp(ii,2);
      (*tt++).set(v,i0,i1,i2,vi.lab);
      numTht[it] = jt;
      it++;      
    }

    for (int ii=0;ii<Thtmp.neb;ii++){
      const Mesh::BorderElement &vi(Thtmp.be(ii));  //const BoundaryEdge &vi(Thtmp(ii));  
      int i0 = nvtmp +  Thtmp.operator()(vi[0]);
      int i1 = nvtmp +  Thtmp.operator()(vi[1]);
      (*bb++).set(v,i0,i1,vi.lab);
      ibe++;
    }
    jt++;
  }
  assert( it==nt ); assert(iv==nv); assert(ibe=nbe);
  if(verbosity) cout << "Popen medit : vertex "<< nv << " triangle "<< nt << " edge " << nbe << endl;  

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

  // 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;
 
  if(nbsol > 0){
    
    vxx.init(datasize);
    vyx.init(datasize);
    vyy.init(datasize);
    
    if( nargs[2] ){
      boolsave= true;
      saveff = GetAny<string *>( (*nargs[2])(stack) );
      int ddim = 2;
    
      for (size_t i=0;i<offset;i++){
	solnbfloat = solnbfloat + l[i].nbfloat;
	//else if( TypTab[i] == 2) solnbfloat = solnbfloat+ddim;
	//else if( TypTab[i] == 3) solnbfloat = solnbfloat+ddim*(ddim+1)/2;
      }
      
      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  Mesh::Vertex & P = Th.vertices[k];
	fx=P.x; fy=P.y;
	fwrite( (unsigned char *) &fx, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &fy, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &(P.lab), WrdSiz,1,popenstream );
	//fprintf(popenstream,"%f %f %i\n",fx,fy,P.lab);
      }
    
      // triangles
      KwdCod = GmfTriangles;
      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 Mesh::Triangle & K(Th.t(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);
      }
      

      // Edges
      KwdCod = GmfEdges;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &nbe, WrdSiz,1,popenstream );
      //fprintf(popenstream,"Edges\n");
      //fprintf(popenstream,"%i\n",nbe);
      for (int k=0; k<nbe; k++) {
	const Mesh::BorderElement & K(Th.be(k));
	int i0=Th.operator()(K[0])+1;
	int i1=Th.operator()(K[1])+1;
	int lab=K.lab;
	fwrite( (unsigned char *) &i0, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &i1, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &lab, WrdSiz,1,popenstream );
	//fprintf(popenstream,"%i %i %i\n",i0,i1,lab);
      }

      // End
      KwdCod = GmfEnd;
      fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
      fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
      //fprintf(popenstream,"End\n");
    }
    else{
      // determination of number solutions associated with a mesh
      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  Mesh::Vertex & P = Th.vertices[k];
	fx=P.x; fy=P.y;
	fprintf(popenstream,"%f %f %i\n",fx,fy,P.lab);
      }
    
      fprintf(popenstream,"Triangles\n");
      fprintf(popenstream,"%i\n",nt);
      for (int k=0; k<nt; k++) {
	const Mesh::Triangle & K(Th.t(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,"Edges\n");
      fprintf(popenstream,"%i\n",nbe);
      for (int k=0; k<nbe; k++) {
	const Mesh::BorderElement & K(Th.be(k));
	int i0=Th.operator()(K[0])+1;
	int i1=Th.operator()(K[1])+1;
	int lab=K.lab;
	fprintf(popenstream,"%i %i %i\n",i0,i1,lab);
      }
      fprintf(popenstream,"End\n");
    }


    // solution with a mesh
    if( nbsol > 0){
      
      if(filebin){
	int cod = 1;
	int NulPos = 0;
	int KwdCod;
	int codtypjm = 1;
	// 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 );
	
	if(order==0){
	  KwdCod = GmfSolAtTriangles;
	  fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &nt, WrdSiz,1,popenstream );

	  printf("SolAtTriangles nt=%i\n",nt);
	}
	if(order==1){
	  KwdCod = GmfSolAtVertices;
	  fwrite( (unsigned char *) &KwdCod, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &NulPos, WrdSiz,1,popenstream );
	  fwrite( (unsigned char *) &nv, WrdSiz,1,popenstream );
	  
	  printf("SolAtVertices nv=%i\n",nv);

	}
	typsol = TypTab[jojo];

	fwrite( (unsigned char *) &codtypjm, WrdSiz,1,popenstream );
	fwrite( (unsigned char *) &typsol, WrdSiz,1,popenstream );
	//printf(popenstream,"%i %i\n",codtypjm,typsol);
      }
      else{
	fprintf(popenstream,"MeshVersionFormatted %i\n",ver);
	fprintf(popenstream,"Dimension %i\n",dimp);
	if(order==0){
	  fprintf(popenstream,"SolAtTriangles\n");
	  fprintf(popenstream,"%i\n",nt);
	}
	if(order==1){
	  fprintf(popenstream,"SolAtVertices\n");
	  fprintf(popenstream,"%i\n",nv);
	}