pos_element.c

cfd求解器使用与gmsh网格的求解

      for(j = 0 ; j < C[i]->NbrNodes ; j++) C[i]->NumNodes[j] = -1 ;
      Cut_PostElement(C[i], GE, PE_L, Index, Depth, Skin, DecomposeInSimplex);
    }

  }
  else{

    Get_InitDofOfElement(&E) ;
    E.GeoElement = GE ;
    E.Num = E.GeoElement->Num ;
    E.Type = E.GeoElement->Type ;
    E.Region = E.GeoElement->Region ;
    Get_NodesCoordinatesOfElement(&E) ;

    for(i = 0 ; i < PE->NbrNodes ; i++){
      if( Skin == 0 && PE->Depth == 1 && 
	  ( DecomposeInSimplex == 0 || E.GeoElement->Type == LINE || 
	    E.GeoElement->Type == TRIANGLE || E.GeoElement->Type == TETRAHEDRON ) ){
	PE->x[i] = E.x[i] ; 
	PE->y[i] = E.y[i] ; 
	PE->z[i] = E.z[i] ;
      }
      else{
	Get_BFGeoElement(&E, PE->u[i], PE->v[i], PE->w[i]) ;
	
	PE->x[i] = PE->y[i] = PE->z[i] = 0. ;
	for (j = 0 ; j < E.GeoElement->NbrNodes ; j++) {
	  PE->x[i] += E.x[j] * E.n[j] ;
	  PE->y[i] += E.y[j] * E.n[j] ;
	  PE->z[i] += E.z[j] * E.n[j] ;
	}
      }      

    }
    
    List_Add(PE_L, &PE);
  }
  
  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  F i l l _ P o s t E l e m e n t                                         */
/* ------------------------------------------------------------------------ */

#define POS_CUT_FILL  Cut_PostElement(PE, GE, PE_L, Index, Depth, 0, DecomposeInSimplex)
#define POS_CUT_SKIN  Cut_PostElement(PE, GE, PE_L, Index, Depth, 1, DecomposeInSimplex)

void Fill_PostElement(struct Geo_Element * GE, List_T * PE_L, 
		      int Index, int Depth, int Skin, List_T * EvaluationPoints_L,
		      int DecomposeInSimplex){

  struct PostElement * PE ;
  int Nbr_EP, i_EP;

  GetDP_Begin("Fill_PostElement");

  if(!Depth){

    PE = Create_PostElement(Index, POINT, 1, 0) ; 
    switch(GE->Type){      
    case POINT       : PE->u[0] = 0.   ; PE->v[0] = 0.   ; PE->w[0] = 0.   ; break ;
    case LINE        : PE->u[0] = 0.   ; PE->v[0] = 0.   ; PE->w[0] = 0.   ; break ;
    case TRIANGLE    : PE->u[0] = 1./3.; PE->v[0] = 1./3.; PE->w[0] = 0.   ; break ;
    case QUADRANGLE  : PE->u[0] = 0.   ; PE->v[0] = 0.   ; PE->w[0] = 0.   ; break ;
    case TETRAHEDRON : PE->u[0] = 0.25 ; PE->v[0] = 0.25 ; PE->w[0] = 0.25 ; break ;
    case HEXAHEDRON  : PE->u[0] = 0.   ; PE->v[0] = 0.   ; PE->w[0] = 0.   ; break ;
    case PRISM       : PE->u[0] = 1./3.; PE->v[0] = 1./3.; PE->w[0] = 0.   ; break ;
#if defined(NEW_PYRAMIDS)
    case PYRAMID     : PE->u[0] = 2./5.; PE->v[0] = 2./5.; PE->w[0] = 1./5.; break ;
#else
    case PYRAMID     : PE->u[0] = 0.   ; PE->v[0] = 0.   ; PE->w[0] = 1./3.; break ;
#endif
    }
    POS_CUT_FILL ;     

  }
  else{

    if(!Skin){

      switch(GE->Type){
      
      case POINT :    
	PE = Create_PostElement(Index, POINT, 1, 1) ; /* node 1 */
	PE->NumNodes[0] = GE->NumNodes[0] ;
	PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
	POS_CUT_FILL ; 
	break ;

      case LINE :    
	PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 1 2 */
	PE->NumNodes[0] = GE->NumNodes[0] ;
	PE->NumNodes[1] = GE->NumNodes[1] ;
	PE->u[0] =-1. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
	PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
	POS_CUT_FILL ; 
	break ;
      
      case TRIANGLE :    
	PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 2 3 */
	PE->NumNodes[0] = GE->NumNodes[0] ;
	PE->NumNodes[1] = GE->NumNodes[1] ;
	PE->NumNodes[2] = GE->NumNodes[2] ;
	PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
	PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
	PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
	POS_CUT_FILL ; 
	break ;
	
      case QUADRANGLE :    
	if(DecomposeInSimplex){
	  PE = Create_PostElement(Index, TRIANGLE, 3, 1); /* nodes 1 2 4 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[1] ;
	  PE->NumNodes[2] = GE->NumNodes[3] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
	  PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
	  PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
	  POS_CUT_FILL;

	  PE = Create_PostElement(Index, TRIANGLE, 3, 1); /* nodes 2 3 4 */
	  PE->NumNodes[0] = GE->NumNodes[1] ;
	  PE->NumNodes[1] = GE->NumNodes[2] ;
	  PE->NumNodes[2] = GE->NumNodes[3] ;
	  PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 0. ; 	      
	  PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 0. ; 	      
	  PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 0. ; 	     
	  POS_CUT_FILL;
	}
	else{
	  if (!EvaluationPoints_L) {
	    PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 2 3 4 */
	    PE->NumNodes[0] = GE->NumNodes[0] ;
	    PE->NumNodes[1] = GE->NumNodes[1] ;
	    PE->NumNodes[2] = GE->NumNodes[2] ;
	    PE->NumNodes[3] = GE->NumNodes[3] ;
	    PE->u[0] = -1. ; PE->v[0] = -1. ; PE->w[0] = 0. ;
	    PE->u[1] =  1. ; PE->v[1] = -1. ; PE->w[1] = 0. ;
	    PE->u[2] =  1. ; PE->v[2] =  1. ; PE->w[2] = 0. ;
	    PE->u[3] = -1. ; PE->v[3] =  1. ; PE->w[3] = 0. ;
	  }
	  else { /* Only for Quadrangles now, to be extended... */
	    Nbr_EP = List_Nbr(EvaluationPoints_L)/3;
	    PE = Create_PostElement(Index, QUADRANGLE, Nbr_EP, 1) ;
	    for (i_EP=0 ; i_EP<Nbr_EP ; i_EP++) {
	      List_Read(EvaluationPoints_L, i_EP*3+0, &PE->u[i_EP]);
	      List_Read(EvaluationPoints_L, i_EP*3+1, &PE->v[i_EP]);
	      List_Read(EvaluationPoints_L, i_EP*3+2, &PE->w[i_EP]);
	    }
	  }
	  POS_CUT_FILL ; 
	}
	break ;
	
      case TETRAHEDRON :
	PE = Create_PostElement(Index, TETRAHEDRON, 4, 1) ; /* nodes 1 2 3 4 */
	PE->NumNodes[0] = GE->NumNodes[0] ;
	PE->NumNodes[1] = GE->NumNodes[1] ;
	PE->NumNodes[2] = GE->NumNodes[2] ;
	PE->NumNodes[3] = GE->NumNodes[3] ;
	PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
	PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ; 	      
	PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 0. ; 	     
	PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ; 	     
	POS_CUT_FILL;
	break ;
	
      case HEXAHEDRON :
	if(DecomposeInSimplex){
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 2 3 6 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[1] ;
	  PE->NumNodes[2] = GE->NumNodes[2] ;
	  PE->NumNodes[3] = GE->NumNodes[5] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] =-1. ; 	     
	  PE->u[3] = 1. ; PE->v[3] =-1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 3 6 7 */    
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[2] ;
	  PE->NumNodes[2] = GE->NumNodes[5] ;
	  PE->NumNodes[3] = GE->NumNodes[6] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ; 	     
	  PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 5 6 7 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[4] ;
	  PE->NumNodes[2] = GE->NumNodes[5] ;
	  PE->NumNodes[3] = GE->NumNodes[6] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ; 	     
	  PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	  
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 3 4 7 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[2] ;
	  PE->NumNodes[2] = GE->NumNodes[3] ;
	  PE->NumNodes[3] = GE->NumNodes[6] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ; 	      
	  PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] =-1. ; 	     
	  PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 5 7 8 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[4] ;
	  PE->NumNodes[2] = GE->NumNodes[6] ;
	  PE->NumNodes[3] = GE->NumNodes[7] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ; 	     
	  PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	  
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 4 7 8 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[3] ;
	  PE->NumNodes[2] = GE->NumNodes[6] ;
	  PE->NumNodes[3] = GE->NumNodes[7] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
	  PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	}
	else{
	  PE = Create_PostElement(Index, HEXAHEDRON, 8, 1) ; /* nodes 1 2 3 4 5 6 7 8 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[1] ;
	  PE->NumNodes[2] = GE->NumNodes[2] ;
	  PE->NumNodes[3] = GE->NumNodes[3] ;
	  PE->NumNodes[4] = GE->NumNodes[4] ;
	  PE->NumNodes[5] = GE->NumNodes[5] ;
	  PE->NumNodes[6] = GE->NumNodes[6] ;
	  PE->NumNodes[7] = GE->NumNodes[7] ;
	  PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
	  PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
	  PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] =-1. ; 	     
	  PE->u[4] =-1. ; PE->v[4] =-1. ; PE->w[4] = 1. ;
	  PE->u[5] = 1. ; PE->v[5] =-1. ; PE->w[5] = 1. ; 	      
	  PE->u[6] = 1. ; PE->v[6] = 1. ; PE->w[6] = 1. ;
	  PE->u[7] =-1. ; PE->v[7] = 1. ; PE->w[7] = 1. ; 	     
	  POS_CUT_FILL;
	}
	break ;
	
      case PRISM :
	if(DecomposeInSimplex){
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 2 3 5 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[1] ;
	  PE->NumNodes[2] = GE->NumNodes[2] ;
	  PE->NumNodes[3] = GE->NumNodes[4] ;
	  PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
	  PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] =-1. ; 	     
	  PE->u[3] = 1. ; PE->v[3] = 0. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	  
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 3 5 6 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[2] ;
	  PE->NumNodes[2] = GE->NumNodes[4] ;
	  PE->NumNodes[3] = GE->NumNodes[5] ;
	  PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
	  PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ; 	     
	  PE->u[3] = 0. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	  
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 4 5 6 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[3] ;
	  PE->NumNodes[2] = GE->NumNodes[4] ;
	  PE->NumNodes[3] = GE->NumNodes[5] ;
	  PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
	  PE->u[1] = 0. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
	  PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ; 	     
	  PE->u[3] = 0. ; PE->v[3] = 1. ; PE->w[3] = 1. ; 	     
	  POS_CUT_FILL;
	}
	else{
	  PE = Create_PostElement(Index, PRISM, 6, 1) ; /* nodes 1 2 3 4 5 6 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[1] ;
	  PE->NumNodes[2] = GE->NumNodes[2] ;
	  PE->NumNodes[3] = GE->NumNodes[3] ;
	  PE->NumNodes[4] = GE->NumNodes[4] ;
	  PE->NumNodes[5] = GE->NumNodes[5] ;
	  PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
	  PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] =-1. ; 	      
	  PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
	  PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ; 	     
	  PE->u[4] = 1. ; PE->v[4] = 0. ; PE->w[4] = 1. ;
	  PE->u[5] = 0. ; PE->v[5] = 1. ; PE->w[5] = 1. ; 	      
	  POS_CUT_FILL;
	}
	break ;
	
      case PYRAMID :
	if(DecomposeInSimplex){
	  PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 2 4 5 */
	  PE->NumNodes[0] = GE->NumNodes[0] ;
	  PE->NumNodes[1] = GE->NumNodes[1] ;
	  PE->NumNodes[2] = GE->NumNodes[3] ;
	  PE->NumNodes[3] = GE->NumNodes[4] ;
#if defined(NEW_PYRAMIDS)
	  PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
	  PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ; 	      
	  PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 0. ;