pos_print.c

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

	  e[NbCut].y[0] = Element.y[abs(NumNodes[0])-1] ;
	  e[NbCut].z[0] = Element.z[abs(NumNodes[0])-1] ;
	  e[NbCut].x[1] = Element.x[abs(NumNodes[1])-1] ;
	  e[NbCut].y[1] = Element.y[abs(NumNodes[1])-1] ;
	  e[NbCut].z[1] = Element.z[abs(NumNodes[1])-1] ;
	  d1 = Plane(A,B,C,D,e[NbCut].x[0],e[NbCut].y[0],e[NbCut].z[0]);
	  d2 = Plane(A,B,C,D,e[NbCut].x[1],e[NbCut].y[1],e[NbCut].z[1]);
	  
	  if(d1*d2 <= 0) {
	    if(d1*d2 < 0.) u = -d2/(d1-d2) ;	    
	    else if(d1 == 0.) u = 1. ;
	    else u = 0. ;
	    e[NbCut].xc = u*e[NbCut].x[0] + (1.-u)*e[NbCut].x[1];
	    e[NbCut].yc = u*e[NbCut].y[0] + (1.-u)*e[NbCut].y[1];
	    e[NbCut].zc = u*e[NbCut].z[0] + (1.-u)*e[NbCut].z[1];	  

	    if(NCPQ_P)
	      xyz2uvwInAnElement(&Element, e[NbCut].xc, e[NbCut].yc, e[NbCut].zc, 
				 &e[NbCut].uc, &e[NbCut].vc, &e[NbCut].wc);	  
	    NbCut++;
	  }
	}
	
	if(NbCut > 3){
	  xcg = ycg = zcg = 0.;
	  for(iCut = 0 ; iCut < NbCut ; iCut++){
	    xcg += e[iCut].xc; ycg += e[iCut].yc; zcg += e[iCut].zc;
	  }
	  xcg /= (double)NbCut; ycg /= (double)NbCut; zcg /= (double)NbCut;
	  DIRZ[0] = A; DIRY[0] = xcg-e[0].xc; 
	  DIRZ[1] = B; DIRY[1] = ycg-e[0].yc; 
	  DIRZ[2] = C; DIRY[2] = zcg-e[0].zc;
	  normvec(DIRZ); normvec(DIRY); prodvec(DIRY,DIRZ,DIRX); normvec(DIRX);
	  XCP = xcg*DIRX[0] + ycg*DIRX[1] + zcg*DIRX[2];
	  YCP = xcg*DIRY[0] + ycg*DIRY[1] + zcg*DIRY[2];	
	  qsort(e,NbCut,sizeof(struct CutEdge), fcmp_Angle);
	}
	
	if(NbCut > 2){
	  iCut = 2;
	  while(iCut < NbCut){
	    if(PostSubOperation_P->Depth > 0){
	      PE = Create_PostElement(iGeo, TRIANGLE, 3, 1) ;
	      PE->x[0] = e[0].xc; PE->x[1] = e[iCut-1].xc; PE->x[2] = e[iCut].xc;
	      PE->y[0] = e[0].yc; PE->y[1] = e[iCut-1].yc; PE->y[2] = e[iCut].yc;
	      PE->z[0] = e[0].zc; PE->z[1] = e[iCut-1].zc; PE->z[2] = e[iCut].zc;
	      PE->u[0] = e[0].uc; PE->u[1] = e[iCut-1].uc; PE->u[2] = e[iCut].uc;
	      PE->v[0] = e[0].vc; PE->v[1] = e[iCut-1].vc; PE->v[2] = e[iCut].vc;
	      PE->w[0] = e[0].wc; PE->w[1] = e[iCut-1].wc; PE->w[2] = e[iCut].wc;
	      LETS_PRINT_THE_RESULT ;
	    }
	    else{
	      PE = Create_PostElement(iGeo, POINT, 1, 0) ;
	      PE->x[0] = (e[0].xc + e[iCut-1].xc + e[iCut].xc) / 3. ;
	      PE->y[0] = (e[0].yc + e[iCut-1].yc + e[iCut].yc) / 3. ;
	      PE->z[0] = (e[0].zc + e[iCut-1].zc + e[iCut].zc) / 3. ;
	      PE->u[0] = (e[0].uc + e[iCut-1].uc + e[iCut].uc) / 3. ;
	      PE->v[0] = (e[0].vc + e[iCut-1].vc + e[iCut].vc) / 3. ;
	      PE->w[0] = (e[0].wc + e[iCut-1].wc + e[iCut].wc) / 3. ;
	      LETS_PRINT_THE_RESULT ;
	    }
	    iCut++;
	  }
	}
	
	if(NbCut == 2){
	  if(PostSubOperation_P->Depth > 0){
	    PE = Create_PostElement(iGeo, LINE, 2, 1) ;
	    PE->x[0] = e[0].xc; PE->x[1] = e[1].xc; 
	    PE->y[0] = e[0].yc; PE->y[1] = e[1].yc; 
	    PE->z[0] = e[0].zc; PE->z[1] = e[1].zc; 
	    PE->u[0] = e[0].uc; PE->u[1] = e[1].uc;
	    PE->v[0] = e[0].vc; PE->v[1] = e[1].vc;
	    PE->w[0] = e[0].wc; PE->w[1] = e[1].wc;
	    LETS_PRINT_THE_RESULT ;
	  }
	  else{
	    PE = Create_PostElement(iGeo, POINT, 1, 0) ;
	    PE->x[0] = (e[0].xc + e[1].xc) / 2. ; 
	    PE->y[0] = (e[0].yc + e[1].yc) / 2. ; 
	    PE->z[0] = (e[0].zc + e[1].zc) / 2. ; 
	    PE->u[0] = (e[0].uc + e[1].uc) / 2. ;
	    PE->v[0] = (e[0].vc + e[1].vc) / 2. ;
	    PE->w[0] = (e[0].wc + e[1].wc) / 2. ;
	    LETS_PRINT_THE_RESULT ;
	  }
	}

      }

    }
    Format_PostFooter(PostSubOperation_P, 0);
    break;
    
  default :
    Msg(GERROR, "Unknown operation in Print OnSection");
    break;
  }

  List_Delete(PE_L) ;
  if(CPQ_P) Free(CumulativeValues);
  for(iCut = 0 ; iCut < NBR_MAX_CUT ; iCut++) Free(e[iCut].Value) ;

  GetDP_End ;
}

#undef NBR_MAX_CUT
#undef LETS_PRINT_THE_RESULT

/* ------------------------------------------------------------------------ */
/*  P o s _ P r i n t O n G r i d                                           */
/* ------------------------------------------------------------------------ */


#define LETS_PRINT_THE_RESULT							\
 PE->x[0] = Current.xp = Current.x ;						\
 PE->y[0] = Current.yp = Current.y ;						\
 PE->z[0] = Current.zp = Current.z ;						\
 if(!NCPQ_P){									\
   for (ts = 0 ; ts < NbTimeStep ; ts++){					\
     PE->Value[0] = CumulativeValues[ts] ;					\
     Format_PostElement(PSO_P->Format, PSO_P->Iso, 0,				\
		        Current.Time, ts, NbTimeStep,				\
                        Current.NbrHar, PSO_P->HarmonicToTime,                  \
			Normal, PE,                                             \
			PSO_P->ChangeOfCoordinates,                             \
			PSO_P->ChangeOfValues);                                 \
   }										\
 }										\
 else{										\
   InWhichElement(Current.GeoData->Grid, NULL, &Element, PSO_P->Dimension,	\
                  Current.x, Current.y, Current.z, &u, &v, &w) ;		\
   Current.Region = Element.Region ;						\
   if(Element.Num != NO_ELEMENT)						\
     PE->Index = Geo_GetGeoElementIndex(Element.GeoElement) ;			\
   else										\
     PE->Index = NO_ELEMENT ;							\
   for (ts = 0 ; ts < NbTimeStep ; ts++) {					\
     Pos_InitAllSolutions(PSO_P->TimeStep_L, ts) ;				\
     Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,		\
                 NULL, &Element, u, v, w, &PE->Value[0]);			\
     if(CPQ_P)									\
       Combine_PostQuantity(PSO_P->CombinationType, Order,			\
                            &PE->Value[0], &CumulativeValues[ts]) ;		\
     Format_PostElement(PSO_P->Format, PSO_P->Iso, 0,				\
                        Current.Time, ts, NbTimeStep,				\
                        Current.NbrHar, PSO_P->HarmonicToTime,                  \
			Normal, PE,                                      	\
			PSO_P->ChangeOfCoordinates,                             \
			PSO_P->ChangeOfValues);                                 \
   }										\
 }

#define ARRAY(i,j,k,t)						\
  Array[ (t) * Current.NbrHar * ((int)N[0]+1) * ((int)N[1]+1) +	\
         (k) * ((int)N[0]+1) * ((int)N[1]+1) + 			\
         (j) * ((int)N[0]+1) +					\
         (i) ]

#define LETS_STORE_THE_RESULT								\
 if(!NCPQ_P){										\
   if(CumulativeValues[0].Type != SCALAR)						\
     Msg(GERROR, "Print OnPlane is not designed for non scalar values with Depth > 1");	\
   else											\
     for (ts = 0 ; ts < NbTimeStep ; ts++)						\
       for(k = 0 ; k < Current.NbrHar ; k++)						\
         ARRAY(i1,i2,k,ts) = (float)CumulativeValues[ts].Val[MAX_DIM*k] ;		\
 }											\
 else{											\
   InWhichElement(Current.GeoData->Grid, NULL, &Element, PSO_P->Dimension,		\
                  Current.x, Current.y, Current.z, &u, &v, &w) ;			\
   Current.Region = Element.Region ;							\
   for (ts = 0 ; ts < NbTimeStep ; ts++) {						\
     Pos_InitAllSolutions(PSO_P->TimeStep_L, ts) ;					\
     Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,			\
                      NULL, &Element, u, v, w, &PE->Value[0]);				\
     if(PE->Value[0].Type != SCALAR)							\
       Msg(GERROR, "Print OnPlane is not designed for non scalar values with Depth > 1");\
     if(CPQ_P)										\
       Combine_PostQuantity(PSO_P->CombinationType, Order,				\
                            &PE->Value[0], &CumulativeValues[ts]) ;			\
     for(k = 0 ; k < Current.NbrHar ; k++)						\
       ARRAY(i1,i2,k,ts) = (float)PE->Value[0].Val[MAX_DIM*k] ;				\
   }											\
 }

void  Pos_PrintOnGrid(struct PostQuantity     *NCPQ_P,
		      struct PostQuantity     *CPQ_P,
		      int                      Order,
		      struct DefineQuantity   *DefineQuantity_P0,
		      struct QuantityStorage  *QuantityStorage_P0,
		      struct PostSubOperation *PSO_P) {

  struct Element       Element ;
  struct Value       * CumulativeValues, Value ;
  struct PostElement * PE , * PE2 ;

  int     i1, i2, i3, k, NbTimeStep, ts ;
  float  *Array = NULL ;
  double  u, v, w, Length, Normal[4] = {0., 0., 0., 0.} ;
  double  X[4], Y[4], Z[4], S[4], N[4], tmp[3];

  GetDP_Begin("Pos_PrintOnGrid");

  Get_InitDofOfElement(&Element) ;

  NbTimeStep = Pos_InitTimeSteps(PSO_P);

  if(CPQ_P){
    Cal_PostCumulativeQuantity(NULL, PSO_P->PostQuantitySupport[Order],
			       PSO_P->TimeStep_L,
			       CPQ_P, DefineQuantity_P0,
			       QuantityStorage_P0, &CumulativeValues);
  }

  Init_SearchGrid(&Current.GeoData->Grid) ;

  Format_PostHeader(PSO_P->Format, PSO_P->Iso, 
		    NbTimeStep, PSO_P->HarmonicToTime,
		    PSO_P->CombinationType, Order,
		    NCPQ_P?NCPQ_P->Name:NULL,
		    CPQ_P?CPQ_P->Name:NULL);

  PE = Create_PostElement(0, POINT, 1, 0) ;

  switch(PSO_P->SubType) {

  case PRINT_ONGRID_0D :
    Current.x = PSO_P->Case.OnGrid.x[0] ;
    Current.y = PSO_P->Case.OnGrid.y[0] ;
    Current.z = PSO_P->Case.OnGrid.z[0] ;
    Normal[0] = Normal[1] = Normal[2] = 0.0 ;
    LETS_PRINT_THE_RESULT ;
    break;

  case PRINT_ONGRID_1D :
    X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
    Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ; 
    Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
    N[0] = PSO_P->Case.OnGrid.n[0] ;
    Normal[1] = Normal[2] = 0.0 ;
    Length = sqrt(DSQU(X[1]-X[0]) + DSQU(Y[1]-Y[0]) + DSQU(Z[1]-Z[0])) ;
    for (i1 = 0 ; i1 <= N[0] ; i1++) {
      S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
      Normal[0] = Length * S[0] ;
      Current.x = X[0] + (X[1] - X[0]) * S[0] ;
      Current.y = Y[0] + (Y[1] - Y[0]) * S[0] ;
      Current.z = Z[0] + (Z[1] - Z[0]) * S[0] ;
      LETS_PRINT_THE_RESULT ;
    }
    break;

  case PRINT_ONGRID_2D :
    X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
    Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
    Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
    X[2] = PSO_P->Case.OnGrid.x[2] ; 
    Y[2] = PSO_P->Case.OnGrid.y[2] ; 
    Z[2] = PSO_P->Case.OnGrid.z[2] ; 

    S[0] = X[1]-X[0]; S[1] = Y[1]-Y[0]; S[2] = Z[1]-Z[0];
    N[0] = X[2]-X[0]; N[1] = Y[2]-Y[0]; N[2] = Z[2]-Z[0];
    prodvec(S,N,Normal);
    Length = sqrt(DSQU(Normal[0])+DSQU(Normal[1])+DSQU(Normal[2]));
    if(!Length){
      Msg(WARNING, "Bad plane (null normal)"); 
      GetDP_End ;
    }
    Normal[0]/=Length ; Normal[1]/=Length ; Normal[2]/=Length ;

    N[0] = PSO_P->Case.OnGrid.n[0] ; N[1] = PSO_P->Case.OnGrid.n[1] ;

    if(PSO_P->Depth > 1)
      Array = (float*)
	Malloc(NbTimeStep*Current.NbrHar*(int)((N[0]+1)*(N[1]+1))*sizeof(float)) ;
    
    for (i1 = 0 ; i1 <= N[0] ; i1++) {
      S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
      for (i2 = 0 ; i2 <= N[1] ; i2++) {
	S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
	Current.x = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[1] ;
	Current.y = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[1] ;
	Current.z = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[1] ;	  
	if(PSO_P->Depth > 1){
	  LETS_STORE_THE_RESULT ;
	}
	else{
	  LETS_PRINT_THE_RESULT ;
	}
      }
      if(PSO_P->Depth < 2 && !Flag_BIN) fprintf(PostStream, "\n");
    }

    if(PSO_P->Depth > 1){
      PE2 = Create_PostElement(0, TRIANGLE, 3, 0);
      PE2->Value[0].Type = PE2->Value[1].Type = PE2->Value[2].Type = SCALAR ;
      for (i1 = 0 ; i1 < N[0] ; i1++) {
	S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
	S[2] = (double)(i1+1) / (double)(N[0] ? N[0] : 1) ;
	for (i2 = 0 ; i2 < N[1] ; i2++) {
	  S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
	  S[3] = (double)(i2+1) / (double)(N[1] ? N[1] : 1) ;
	  PE2->x[0] = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[1] ;
	  PE2->y[0] = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[1] ;
	  PE2->z[0] = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[1] ;	  
	  PE2->x[1] = X[0] + (X[1] - X[0]) * S[2] + (X[2] - X[0]) * S[1] ;
	  PE2->y[1] = Y[0] + (Y[1] - Y[0]) * S[2] + (Y[2] - Y[0]) * S[1] ;
	  PE2->z[1] = Z[0] + (Z[1] - Z[0]) * S[2] + (Z[2] - Z[0]) * S[1] ;	  
	  PE2->x[2] = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[3] ;
	  PE2->y[2] = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[3] ;
	  PE2->z[2] = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[3] ;	  
	  for (ts = 0 ; ts < NbTimeStep ; ts++){
	    for(k = 0 ; k < Current.NbrHar ; k++){
	      PE2->Value[0].Val[MAX_DIM*k] = ARRAY(i1,i2,k,ts) ;
	      PE2->Value[1].Val[MAX_DIM*k] = ARRAY(i1+1,i2,k,ts) ;
	      PE2->Value[2].Val[MAX_DIM*k] = ARRAY(i1,i2+1,k,ts) ;
	    }
	    Format_PostElement(PSO_P->Format, PSO_P->Iso, 0,
			       Current.Time, ts, NbTimeStep,
			       Current.NbrHar, PSO_P->HarmonicToTime, 
			       Normal, PE2,
			       PSO_P->ChangeOfCoordinates,
			       PSO_P->ChangeOfValues);
	  }
	  PE2->x[0] = X[0] + (X[1] - X[0]) * S[2] + (X[2] - X[0]) * S[3] ;
	  PE2->y[0] = Y[0] + (Y[1] - Y[0]) * S[2] + (Y[2] - Y[0]) * S[3] ;
	  PE2->z[0] = Z[0] + (Z[1] - Z[0]) * S[2] + (Z[2] - Z[0]) * S[3] ;	  
	  tmp[0] = PE2->x[1]; PE2->x[1] = PE2->x[2]; PE2->x[2] = tmp[0];
	  tmp[1] = PE2->y[1]; PE2->y[1] = PE2->y[2]; PE2->y[2] = tmp[1];
	  tmp[2] = PE2->z[1]; PE2->z[1] = PE2->z[2]; PE2->z[2] = tmp[2];
	  for (ts = 0 ; ts < NbTimeStep ; ts++){
	    for(k = 0 ; k < Current.NbrHar ; k++){
	      PE2->Value[0].Val[MAX_DIM*k] = ARRAY(i1+1,i2+1,k,ts) ;
	      PE2->Value[1].Val[MAX_DIM*k] = ARRAY(i1,i2+1,k,ts) ;
	      PE2->Value[2].Val[MAX_DIM*k] = ARRAY(i1+1,i2,k,ts) ;
	    }
	    Format_PostElement(PSO_P->Format, PSO_P->Iso, 0,
			       Current.Time, ts, NbTimeStep,
			       Current.NbrHar, PSO_P->HarmonicToTime, 
			       Normal, PE2, 
			       PSO_P->ChangeOfCoordinates,
			       PSO_P->ChangeOfValues);
	  }
	}
      }
      Destroy_PostElement(PE2) ;
      Free(Array) ;
    }
    break;

  case PRINT_ONGRID_3D :
    X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ; 
    Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
    Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
    X[2] = PSO_P->Case.OnGrid.x[2] ; X[3] = PSO_P->Case.OnGrid.x[3] ;
    Y[2] = PSO_P->Case.OnGrid.y[2] ; Y[3] = PSO_P->Case.OnGrid.y[3] ;
    Z[2] = PSO_P->Case.OnGrid.z[2] ; Z[3] = PSO_P->Case.OnGrid.z[3] ;
    N[0] = PSO_P->Case.OnGrid.n[0] ; 
    N[1] = PSO_P->Case.OnGrid.n[1] ;
    N[2] = PSO_P->Case.OnGrid.n[2] ;
    Normal[0] = Normal[1] = Normal[2] = 0.0 ;
    for (i1 = 0 ; i1 <= N[0] ; i1++) {
      S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
      for (i2 = 0 ; i2 <= N[1] ; i2++) {
	S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
	for (i3 = 0 ; i3 <= N[2] ; i3++) {
	  S[2] = (double)i3 / (double)(N[2] ? N[2] : 1) ;
	  Current.x = X[0] + (X[1]-X[0])*S[0] + (X[2]-X[0])*S[1] + (X[3]-X[0])*S[2] ;
	  Current.y = Y[0] + (Y[1]-Y[0])*S[0] + (Y[2]-Y[0])*S[1] + (Y[3]-Y[0])*S[2] ;
	  Current.z = Z[0] + (Z[1]-Z[0])*S[0] + (Z[2]-Z[0])*S[1] + (Z[3]-Z[0])*S[2] ;
	  LETS_PRINT_THE_RESULT ;
	}
	if(!Flag_BIN) fprintf(PostStream, "\n");
      }
      if(!Flag_BIN) fprintf(PostStream, "\n\n");
      /*  two blanks lines for -index in gnuplot  */
    }
    break;

  case PRINT_ONGRID_PARAM :
    for (i1 = 0 ; i1 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[0]) ; i1++) {
      List_Read(PSO_P->Case.OnParamGrid.ParameterValue[0], i1, &Current.a) ;
      for (i2 = 0 ; i2 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[1]) ; i2++) {
	List_Read(PSO_P->Case.OnParamGrid.ParameterValue[1], i2, &Current.b) ;
	for (i3 = 0 ; i3 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[2]) ; i3++) {
	  List_Read(PSO_P->Case.OnParamGrid.ParameterValue[2], i3, &Current.c) ;
	  Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[0],
				       NULL, 0., 0., 0., &Value) ; 
	  Current.x = Value.Val[0];
	  Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[1],
				       NULL, 0., 0., 0., &Value) ; 
	  Current.y = Value.Val[0];
	  Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[2],