mesh3dn.cpp

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

    
    for (int k=0; k<nbe; k++) {
      int i[4],lab;
      BorderElement & K(this->borderelements[k]);
      f >> i[0] >> i[1] >> i[2]  >> lab;
      K.set(this->vertices,i,lab);
      mesb += K.mesure();	    
      
    }
    f >> s;
    ffassert( s== Gsend);
  }
  
  
  
  int Mesh3::Save(const string & filename) const
  {
    int ver = GmfFloat, outm;
    if ( !(outm = GmfOpenMesh(filename.c_str(),GmfWrite,ver,3)) ) {
      cerr <<"  -- Mesh3::Save  UNABLE TO OPEN  :"<< filename << endl;
      return(1);
    }
    
    float fx,fy,fz;
    GmfSetKwd(outm,GmfVertices,this->nv);
    for (int k=0; k<nv; k++) {
      const  Vertex & P = this->vertices[k];
      GmfSetLin(outm,GmfVertices,fx=P.x,fy=P.y,fz=P.z,P.lab);
    }
    
    if(nt != 0){
      GmfSetKwd(outm,GmfTetrahedra,nt);
      for (int k=0; k<nt; k++) {
	const Element & K(this->elements[k]);
	int i0=this->operator()(K[0])+1;
	int i1=this->operator()(K[1])+1;
	int i2=this->operator()(K[2])+1;
	int i3=this->operator()(K[3])+1;
	int lab=K.lab;
	GmfSetLin(outm,GmfTetrahedra,i0,i1,i2,i3,lab);
      }
    }
    
    GmfSetKwd(outm,GmfTriangles,nbe);
    for (int k=0; k<nbe; k++) {
      const BorderElement & K(this->borderelements[k]);
      int i0=this->operator()(K[0])+1;
      int i1=this->operator()(K[1])+1;
      int i2=this->operator()(K[2])+1;
      int lab=K.lab;
      GmfSetLin(outm,GmfTriangles,i0,i1,i2,lab);
    }
    
    GmfCloseMesh(outm);
    return (0);
    
  }

   int Mesh3::SaveSurface(const string & filename) const
  {
    int ver = GmfFloat, outm;
    if ( !(outm = GmfOpenMesh(filename.c_str(),GmfWrite,ver,3)) ) {
      cerr <<"  -- Mesh3::Save  UNABLE TO OPEN  :"<< filename << endl;
      return(1);
    }

    // Number of Vertex in the surface
    int *v_num_surf=new int[nv];
    int *liste_v_num_surf=new int[nv];
    for (int k=0; k<nv; k++){ 
      v_num_surf[k]=-1;
      liste_v_num_surf[k]=0;
    }
    // Search Vertex on the surface
    int nbv_surf=0;
    for (int k=0; k<nbe; k++) {
      const BorderElement & K(this->borderelements[k]);     
      for(int jj=0; jj<3; jj++){
	int i0=this->operator()(K[jj]);
	if( v_num_surf[i0] == -1 ){
	  v_num_surf[i0] = nbv_surf;
	  liste_v_num_surf[nbv_surf]= i0;
	  nbv_surf++;
	}
      }
    }

    float fx,fy,fz;
    GmfSetKwd(outm,GmfVertices,nbv_surf);
    for (int k=0; k<nbv_surf; k++) {
      int k0 = liste_v_num_surf[k];
      const  Vertex & P = this->vertices[k0];
      GmfSetLin(outm,GmfVertices,fx=P.x,fy=P.y,fz=P.z,P.lab);
    }
    
    GmfSetKwd(outm,GmfTriangles,nbe);
    for (int k=0; k<nbe; k++) {
      const BorderElement & K(this->borderelements[k]);
      int i0=v_num_surf[this->operator()(K[0])]+1;
      int i1=v_num_surf[this->operator()(K[1])]+1;
      int i2=v_num_surf[this->operator()(K[2])]+1;
      int lab=K.lab;

      assert( i0-1 < nbv_surf &&  i1-1 < nbv_surf &&  i2-1 < nbv_surf );
      assert( 0<i0 && 0<i1 && 0<i2 );

      GmfSetLin(outm,GmfTriangles,i0,i1,i2,lab);
    }
    
    GmfCloseMesh(outm);

    delete [ ] v_num_surf;
    delete [ ] liste_v_num_surf;

    return (0); 
  }


  int Mesh3::SaveSurface(const string & filename1,const string & filename2) const
  {    
    // Number of Vertex in the surface
    int *v_num_surf=new int[nv];
    int *liste_v_num_surf=new int[nv];
    for (int k=0; k<nv; k++){ 
      v_num_surf[k]=-1;
      liste_v_num_surf[k]=0;
    }
    // Search Vertex on the surface
    int nbv_surf=0;
    for (int k=0; k<nbe; k++) {
      const BorderElement & K(this->borderelements[k]);     
      for(int jj=0; jj<3; jj++){
	int i0=this->operator()(K[jj]);
	if( v_num_surf[i0] == -1){
	  v_num_surf[i0] = nbv_surf;
	  nbv_surf++;
	}
      }
    }

    // file .points
    FILE *fpoints = fopen(filename1.c_str(),"w");
    fprintf(fpoints,"%i\n",nbv_surf);
    
    for (int k=0; k<nbv_surf; k++) {
      int k0 = liste_v_num_surf[k];
      const  Vertex & P = this->vertices[k];
      fprintf(fpoints,"%f %f %f %i\n",P.x,P.y,P.z,P.lab);
    }
    fclose(fpoints);
    
    // file .faces
    FILE *ffaces = fopen(filename2.c_str(),"w");
    fprintf(ffaces,"%i\n",nbe);
    for (int k=0; k<nbe; k++) {
      const BorderElement & K(this->borderelements[k]);
      int i0=this->operator()(K[0]);
      int i1=this->operator()(K[1]);
      int i2=this->operator()(K[2]);
      int lab=K.lab;
      int label0= this->vertices[i0].lab; 
      int label1= this->vertices[i1].lab; 
      int label2= this->vertices[i2].lab;
      //GmfSetLin(outm,GmfTriangles,i0,i1,i2,lab);
      int nature=3;
      int i0v=v_num_surf[i0]+1;
      int i1v=v_num_surf[i1]+1;
      int i2v=v_num_surf[i2]+1;
      assert( i0v-1 < nbv_surf &&  i1v-1 < nbv_surf &&  i2v-1 < nbv_surf );
      assert( 0<i0v && 0<i1v && 0<i2v );

      fprintf(ffaces,"%i %i %i %i %i %i %i %i\n", nature, i0v, i1v, i2v, lab, label0, label1, label2);
    }
    fclose(ffaces);
    
    delete [ ] v_num_surf;
    delete [ ] liste_v_num_surf;
    return (0);  
  }


  Mesh3::Mesh3(int nnv, int nnt, int nnbe, Vertex3 *vv, Tet *tt, Triangle3 *bb)
  {
  
    nv = nnv;
    nt = nnt;
    nbe =nnbe;
    
    vertices = vv;
    elements = tt;
    borderelements = bb;
    
    mes=0.;
    mesb=0.;
    
    for (int i=0;i<nt;i++)  
      mes += this->elements[i].mesure();
    
    for (int i=0;i<nbe;i++)  
      mesb += this->be(i).mesure();  
    
    
    //if(nnt !=0){		      
    //cout << "action  sur le maillage" << endl;
    //BuildBound();
    //BuildAdj();
    //Buildbnormalv();  
    //BuildjElementConteningVertex();
    //BuildGTree();
    //decrement();    
    //}
    
    if(verbosity>1)
      cout << "  -- End of read: mesure = " << mes << " border mesure " << mesb << endl;  
    
  }
  
  Mesh3::Mesh3(int nnv, int nnbe, Vertex3 *vv, Triangle3 *bb)
  {
  
    nv = nnv;
    nbe =nnbe;
    
    vertices = vv;
    borderelements = bb;
    
    mes=0.;
    mesb=0.;
    
    for (int i=0;i<nbe;i++)  
      mesb += this->be(i).mesure();  
    
    if(verbosity>1)
      cout << "  -- End of read: mesure = " << mes << " border mesure " << mesb << endl;  
  }

  void Mesh3::flipSurfaceMesh3(int surface_orientation)
  {
    /* inverse the orientation of the surface if necessary*/
    /* and control that all surfaces are oriented in the same way*/
    int nbflip=0;
    for (int i=0;i<this->nbe;i++)
      { 
	double mes_triangle3= this->be(i).mesure();
	
	if( surface_orientation*mes_triangle3 < 0){
	  const Triangle3 &K( this->be(i) );
	  int iv[3];       
	  
	  iv[0] = this->operator()(K[0]);
	  iv[1] = this->operator()(K[1]);
	  iv[2] = this->operator()(K[2]);
	  
	  int iv_temp=iv[1];
	  iv[1]=iv[2];
	  iv[2]=iv_temp;
	  this->be(i).set( this->vertices, iv, K.lab ) ;
	  nbflip++;
	}
      }
    assert(nbflip==0 || nbflip== this->nbe); 
  }


  int  signe_permutation(int i0,int i1,int i2,int i3)
  {
    int p=1;
    if(i0>i1) Exchange(i0,i1), p = -p;
    if(i0>i2) Exchange(i0,i2), p = -p;
    if(i0>i3) Exchange(i0,i3), p = -p;
    if(i1>i2) Exchange(i1,i2), p = -p;
    if(i1>i3) Exchange(i1,i3), p = -p;
    if(i2>i3) Exchange(i2,i3), p = -p;
    return p;
  }


  int  WalkInTet(const Mesh3 & Th,int it, R3 & Phat,const R3 & U, R & dt)
  {
      bool ddd=verbosity>200;
    R lambda[4];
    Phat.toBary(lambda);
    if(ddd) cout << "\n\n\n   WT: "  << Phat << " :  "  << lambda[0] << " " <<lambda[1] <<" " <<lambda[2] << " " <<lambda[3] << " u = "<< U << " dt " << dt <<endl;
#ifndef NDEBUG      
      for(int i=0;i<4;++i)
      assert(lambda[i]<1.000001 && lambda[i]>-0.0000001);
#endif
    typedef R3 Rd;
    const Mesh3::Element & T(Th[it]);
    const int nve = 4;
    const Rd  Q[nve]={(const R3) T[0],(const R3) T[1],(const R3) T[2],(const R3) T[3]};
    
    Rd P  =T(Phat);
    
    //  cout << " " << u << " " << v ;
    Rd PF = P + U*dt;
    
    //  couleur(15);MoveTo( P); LineTo( PF);
    R l[nve];
    double Det=T.mesure()*6;
    l[0] = det(PF  ,Q[1],Q[2],Q[3]);
    l[1] = det(Q[0],PF  ,Q[2],Q[3]); 
    l[2] = det(Q[0],Q[1],PF  ,Q[3]); 
    l[3] = Det - l[0]-l[1]-l[2];
    l[0] /= Det;
    l[1] /= Det;
    l[2] /= Det;
    l[3] /= Det;
     if(ddd)  cout << "\t\t\tWT " << it << ", " << Phat << ",  " << PF
		   << " :  "  << l[0] << " " <<l[1] <<" " <<l[2] << " " <<l[3] 
	           << " == " << det(Q[0],Q[1],Q[2],PF)/Det
		   << " :  "  << lambda[0] << " " <<lambda[1] <<" " <<lambda[2] << " " <<lambda[3] 
	             
		   <<endl ;
		  
    const R eps = 1e-8;
    int neg[nve],k=0;
    int kk=-1;
    if (l[0]>-eps && l[1]>-eps && l[2]>-eps && l[3]>-eps) 
      {
	dt =0;
	Phat=R3(l+1);

      }
    else 
      {
	// on regarde de les reelement negatif 
        // on ne veut pas des points sur les faces.
        // car sinon il va y avoir un probleme ans on va projete sur la face
	//  et remettre le point dans le tetraedre.
	if (l[0]<=-eps ) neg[k++]=0;
	if (l[1]<=-eps ) neg[k++]=1;
	if (l[2]<=-eps ) neg[k++]=2;
	if (l[3]<=-eps ) neg[k++]=3;
	
	R eps1 = T.mesure()   * 1.e-5;
	   if(ddd)  cout << " k= " << k << endl;
	if (k==3) //  3 face de sortie possible 
	  {
	    // let j be the vertex beetween the 3 faces 
	    int j = 6-neg[0]-neg[1]-neg[2]; //  sommet intersection des 3 faces.
	    int i0 = Tet::nvface[j][0];
	    int i1 = Tet::nvface[j][1];
	    int i2 = Tet::nvface[j][2];
	     if(ddd)  cout << "  -------- " << j << " " << i0 << " " << i1 << " " << i2  << endl;
	    //  le tet i0,i1,i2,j est positif. 
	    assert(signe_permutation(i0,i1,i2,j)==1);
	    // 
	    R v0= det(Q[i0],Q[j],P,PF); 
	    R v1= det(Q[i1],Q[j],P,PF); 
	    R v2= det(Q[i2],Q[j],P,PF); 
	     if(ddd)   cout << "\t\t\t " << j << " v0123 =" << v0 << " "<< v1 << " " << v2 << endl;
	    if( v0 > eps && v1 < -eps ) 
	      kk= i1 ;// on sort par la face j i0, j1
	    else if( v1 > eps && v2 < -eps ) 
	      kk= i0 ;
	    else if( v2 > eps && v0 < -eps ) 
	      kk= i1 ;
	    else 
	      {  // on ne sort pas par une face 
		int nul[3], nn=0, mm=3;
		  if (Abs(v0) <=eps) nul[nn++]=i0; else nul[--mm]=i0;
		  if (Abs(v1) <=eps) nul[nn++]=i1; else nul[--mm]=i1;
		  if (Abs(v2) <=eps) nul[nn++]=i2; else nul[--mm]=i2;
		assert(nn>0); 
		if(nn == 1) // on sort par l'arete nul[0] entre le face   nul[1] et nul[2]
		  kk =  nul[1+(rand()/(RAND_MAX/2))%2];		  
		else // on sort par le sommet j.  on choisi la face alleatoirement 
		  kk = nul[(rand()/(RAND_MAX/3))%3];
		  
	      }
	  }
	else if (k==2)
	  {
	    //  numero des l'arete entre les 2 faces
	    int i0=neg[0],i1=neg[1];
	    int e = i0 + i1 - (i0==0); 
	    // on a:
	    //   e      =     0        1       2      3        4      5
	    //  (i0,i1) =   (0,1)  , (0,2), (0,3) , (1,2)  , (1,3), (2,3) 
	    // avec i0,i1 sont les sommets qui ne sont pas dans l'arete
	    int   jj0[6] = {2,3,1,0,2,0};
	    int   jj1[6] = {3,1,2,3,0,1};
	    int j0 = jj0[e];
	    int j1 = jj1[e];
	     if(ddd)   cout << " e " << e << " i0 " << i0 << " " << i1 << " j0 =" << j0 << " " << j1 << endl;
	    // le tet  j0,j1,i0,i1  doit est positif (ie. la pemutation est positive)
	    // de meme  i0,i1,j0,j1
	    assert(signe_permutation(j0,j1,i0,i1)==1);
	    R v0= det(Q[j0],Q[j1],P,PF); 
            if(ddd) cout << " v0 =" << v0 <<endl;
	    if( Abs(v0) < eps  ) 
	      {
	      // on sort par l'arete  j0,j1
	      // on choisi aleatoirement la face de sortie 
	      kk = (rand()/(RAND_MAX/2)) ? i0 : i1; 
	      if(ddd)
		  cout << " rand choose  2 :  " << kk << endl;
	      }
	    else 
	      kk= v0 >0 ? i0 : i1; // Attention dyslexie ici durdur FH....
	    
	  }
	else if (k==1) //  une face possible de sortie (cas simple)
	  kk = neg[0];
	
	if(kk>=0)
	  {
	    R d=lambda[kk]-l[kk];
	    if ( l[kk] )
	     {
	    throwassert(d);
	    R coef =  lambda[kk]/d;
	    R coef1 = 1-coef;
	    dt        = dt*coef1;
	    lambda[0] = lambda[0]*coef1 + coef *l[0];
	    lambda[1] = lambda[1]*coef1 + coef *l[1];
	    lambda[2] = lambda[2]*coef1 + coef *l[2];
	    lambda[3] = lambda[3]*coef1 + coef *l[3];
            if(ddd) cout << "   \t\t -> kk=" << kk << " d=" << d << " , l= "<< lambda[0]  << " " 
			 <<lambda[1] << " " <<lambda[2] << " " << lambda[3] << endl;
	    lambda[kk] =0;
	     }
	    else // on ne bouge pas on resort 
	      {
	       if(ddd)  cout << "            WT : on ne bouge pas on resort \n";
	      return kk;
	      }
	      
	  }
      }
    //  on remet le point dans le tet. 
    int jj=0;
    R lmx=lambda[0];
    if (lmx<lambda[1])  jj=1, lmx=lambda[1];
    if (lmx<lambda[2])  jj=2, lmx=lambda[2];
    if (lmx<lambda[3])  jj=3, lmx=lambda[3];
    if(lambda[0]<0) lambda[jj] += lambda[0],lambda[0]=0;
    if(lambda[1]<0) lambda[jj] += lambda[1],lambda[1]=0;
    if(lambda[2]<0) lambda[jj] += lambda[2],lambda[2]=0;
    if(lambda[3]<0) lambda[jj] += lambda[3],lambda[3]=0;
    Phat=R3(lambda+1);
    if(ddd) cout  << "\t\t\t -> "<< dt << " : "  << Phat << ", " << kk << " jj= "<< jj << " "<< lmx << endl; 
    assert(kk<0 || lambda[kk]==0);
    return kk;
  }        
	      
    
} //   End  namespace Fem2D