pos_print.c

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

	      Current.z = PE->z[iNode] ;
	      Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,
			       NULL, &Element,
			       PE->u[iNode], PE->v[iNode], PE->w[iNode], &PE->Value[iNode]);
	    if(CPQ_P)
	      Combine_PostQuantity(PostSubOperation_P->CombinationType, Order,
				   &PE->Value[iNode], &CumulativeValues[iTime]) ;
	    }      

	    if(!Store)
	      Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso, 0,
				 Current.Time, iTime, NbrTimeStep,
				 Current.NbrHar, PostSubOperation_P->HarmonicToTime,
				 NULL, PE, 
				 PostSubOperation_P->ChangeOfCoordinates,
				 PostSubOperation_P->ChangeOfValues);
	  }
	}
      }
      
      if(!Store) Destroy_PostElement(PE) ;

    }
  } /* for iGeo */

  /* Perform Smoothing */
  
  if(PostSubOperation_P->Smoothing && !InteractiveInterrupt){

    Msg(INFO, "Smoothing (phase 1)");

    xyzv_T = Tree_Create(sizeof(struct xyzv), fcmp_xyzv);
    
    for (iPost = 0 ; iPost < NbrPost ; iPost++){ 
      PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
      for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++) {
	xyzv.x = PE->x[iNode];
	xyzv.y = PE->y[iNode];
	xyzv.z = PE->z[iNode];
	if((xyzv_P = (struct xyzv*)Tree_PQuery(xyzv_T, &xyzv))){
	  x1 = (double)(xyzv_P->nboccurences)/ (double)(xyzv_P->nboccurences + 1.);
	  x2 = 1./(double)(xyzv_P->nboccurences + 1);
	  Cal_AddMultValue2(&xyzv_P->v, x1, &PE->Value[iNode], x2);
	  xyzv_P->nboccurences++;
	}
	else{
	  Cal_CopyValue(&PE->Value[iNode],&xyzv.v);
	  xyzv.nboccurences = 1;
	  Tree_Add(xyzv_T, &xyzv);
	}
      }
    }

    Msg(INFO, "Smoothing (phase 2)");

    for(iPost = 0 ; iPost < NbrPost ; iPost++){ 
      PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
      for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++) {
	xyzv.x = PE->x[iNode];
	xyzv.y = PE->y[iNode];
	xyzv.z = PE->z[iNode];
	if((xyzv_P = (struct xyzv*)Tree_PQuery(xyzv_T, &xyzv))){
	  Cal_CopyValue(&xyzv_P->v, &PE->Value[iNode]);
	}
	else
	  Msg(WARNING, "Node (%g,%g,%g) not found", xyzv.x, xyzv.y, xyzv.z);
      }
    }
    
    Tree_Delete(xyzv_T);

  } /* if Smoothing */

  /* Perform Adaption */

  if(PostSubOperation_P->Adapt && !InteractiveInterrupt){

    if(!Current.GeoData->H)
      Current.GeoData->H = (double*)Malloc((NbrGeo+2)*sizeof(double)) ;

    if(!Current.GeoData->P)
      Current.GeoData->P = (double*)Malloc((NbrGeo+2)*sizeof(double)) ;
    
    Error = (double*)Malloc((NbrGeo+1)*sizeof(double)) ;

    /* All elements are perfect... */
    for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++){
      Element.GeoElement = Geo_GetGeoElement(iGeo) ;
      Element.Num = Element.GeoElement->Num ;
      Element.Type = Element.GeoElement->Type ;
      Element.Region = Element.GeoElement->Region ;
      Get_NodesCoordinatesOfElement(&Element) ;
    
      Current.GeoData->H[iGeo+1] = Cal_MaxEdgeLength(&Element) ;
      Current.GeoData->P[iGeo+1] = 1. ;
      Error[iGeo+1] = PostSubOperation_P->Target ;
    }

    /* ...except those we want to optimize */
    for(iPost = 0 ; iPost < NbrPost ; iPost++){ 
      PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
      Error[PE->Index+1] = 0. ;
      for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++)
	Error[PE->Index+1] += PE->Value[iNode].Val[0] ;
      Error[PE->Index+1] /= (double)PE->NbrNodes ;
    }

    Adapt (NbrGeo, PostSubOperation_P->Adapt, 
	   PostSubOperation_P->Dimension, Error,
	   Current.GeoData->H, Current.GeoData->P, 
	   PostSubOperation_P->Target);

    /* Clean up the interpolation orders to fit to what's available */
    if(List_Nbr(PostSubOperation_P->Value_L)){
      for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++){
	for(jj = List_Nbr(PostSubOperation_P->Value_L)-1 ; jj >= 0  ; jj--){
	  d = *(double*)List_Pointer(PostSubOperation_P->Value_L, jj);
	  if(Current.GeoData->P[iGeo+1] > d || jj == 0){
	    Current.GeoData->P[iGeo+1] = d ;
	    break ;
	  }
	}
      }
    }
  } /* if Adapt */

  /* Print everything if we are in Store mode */

  if(Store && !InteractiveInterrupt){

    /* Sort the elements */
    
    switch(PostSubOperation_P->Sort){
    case SORT_BY_POSITION : List_Sort(PostElement_L, fcmp_PostElement) ; break ;
    case SORT_BY_CONNECTIVITY : Sort_PostElement_Connectivity(PostElement_L) ; break ;
    }

    Dummy[0] = Dummy[1] = Dummy[2] = Dummy[3] = Dummy[4] = 0. ;

    for(iPost = 0 ; iPost < NbrPost ; iPost++){ 
      PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);

      /* Get the values from adaption */
      if(PostSubOperation_P->Adapt){
	Element.GeoElement = Geo_GetGeoElement(PE->Index) ;

	Dummy[0] = Element.GeoElement->Num ;
	Dummy[1] = Error[PE->Index+1] ;
	Dummy[2] = Current.GeoData->H[PE->Index+1] ;
	Dummy[3] = Current.GeoData->P[PE->Index+1] ;
	Dummy[4] = iPost ? 0 : NbrPost ;
	
	for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
	  PE->Value[iNode].Type = SCALAR ;
	  if(PostSubOperation_P->Adapt == H1 ||
	     PostSubOperation_P->Adapt == H2)
	    PE->Value[iNode].Val[0] = Dummy[2] ;
	  else
	    PE->Value[iNode].Val[0] = Dummy[3] ;
	}
      }

      /* Compute curvilinear coord if connection sort */
      if(PostSubOperation_P->Sort == SORT_BY_CONNECTIVITY){
	Dummy[0] = Dummy[1] ;
	Dummy[1] = Dummy[0] + sqrt(DSQU(PE->x[1]-PE->x[0])+
				   DSQU(PE->y[1]-PE->y[0])+
				   DSQU(PE->z[1]-PE->z[0])) ;
	Dummy[2] = PE->v[0] ;
	Dummy[3] = -1. ;
      }

      Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso, 1,
			 Current.Time, 0, 1, 
			 Current.NbrHar, PostSubOperation_P->HarmonicToTime,
			 Dummy, PE, 
			 PostSubOperation_P->ChangeOfCoordinates,
			 PostSubOperation_P->ChangeOfValues);
    }
  }

  Format_PostFooter(PostSubOperation_P, Store);

  if(Store)
    for(iPost = 0 ; iPost < NbrPost ; iPost++){ 
      PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
      Destroy_PostElement(PE) ;
    }

  List_Delete(PostElement_L);

  if(CPQ_P) Free(CumulativeValues);
  if(PostSubOperation_P->Adapt) Free(Error) ;

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  P o s _ P r i n t O n S e c t i o n                                     */
/* ------------------------------------------------------------------------ */

double Plane(double a, double b, double c, double d, 
	     double x, double y, double z){

  GetDP_Begin("Plane");

  GetDP_Return(a*x+b*y+c*z+d);
}

static double DIRX[3], DIRY[3], DIRZ[3], XCP, YCP ;

int fcmp_Angle (const void *a, const void *b){
  struct CutEdge *q , *w;
  double x1,y1,x2,y2,ang1,ang2;

  q = (struct CutEdge*)a;
  w = (struct CutEdge*)b;
  
  x1 = q->xc*DIRX[0] + q->yc*DIRX[1] + q->zc*DIRX[2];
  y1 = q->xc*DIRY[0] + q->yc*DIRY[1] + q->zc*DIRY[2];
  x2 = w->xc*DIRX[0] + w->yc*DIRX[1] + w->zc*DIRX[2];
  y2 = w->xc*DIRY[0] + w->yc*DIRY[1] + w->zc*DIRY[2];
  
  ang1 = atan2(y1-YCP,x1-XCP);
  ang2 = atan2(y2-YCP,x2-XCP);

  if(ang1>ang2)return 1;
  return -1;
}

void prodvec (double *a , double *b , double *c){

  GetDP_Begin("prodvec");

  c[0] = a[1]*b[2]-a[2]*b[1];
  c[1] = a[2]*b[0]-a[0]*b[2];
  c[2] = a[0]*b[1]-a[1]*b[0];

  GetDP_End ;
}

void normvec(double *a){
  double mod;

  GetDP_Begin("normvec");

  mod = sqrt(SQU(a[0])+SQU(a[1])+SQU(a[2]));
  a[0]/=mod;
  a[1]/=mod;
  a[2]/=mod;

  GetDP_End ;
}

#define NBR_MAX_CUT 10

#define LETS_PRINT_THE_RESULT							\
  List_Reset(PE_L);								\
  if(PostSubOperation_P->Depth < 2)						\
    List_Add(PE_L, &PE) ;							\
  else										\
    Cut_PostElement(PE, Element.GeoElement, PE_L, PE->Index,			\
		    PostSubOperation_P->Depth, 0, 1) ;				\
  for(iPost = 0 ; iPost < List_Nbr(PE_L) ; iPost++){				\
    PE = *(struct PostElement **)List_Pointer(PE_L, iPost) ;			\
    for(iTime = 0 ; iTime < NbTimeStep ; iTime++){				\
      Pos_InitAllSolutions(PostSubOperation_P->TimeStep_L, iTime) ;     	\
      for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){				\
	if(NCPQ_P){								\
	  Current.x = PE->x[iNode] ;						\
	  Current.y = PE->y[iNode] ;						\
	  Current.z = PE->z[iNode] ;						\
	  Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,	\
		      NULL, &Element, PE->u[iNode], PE->v[iNode], PE->w[iNode],	\
		      &PE->Value[iNode]);					\
	  if(CPQ_P)								\
             Combine_PostQuantity(PostSubOperation_P->CombinationType, Order,	\
				  &PE->Value[iNode], &CumulativeValues[iTime]) ;\
	}									\
	else									\
	  Cal_CopyValue(&CumulativeValues[iTime],&PE->Value[iNode]);		\
      }										\
      Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso,0, \
			 Current.Time, iTime, NbTimeStep,			\
			 Current.NbrHar, PostSubOperation_P->HarmonicToTime,	\
			 NULL, PE,                                              \
			 PostSubOperation_P->ChangeOfCoordinates,               \
			 PostSubOperation_P->ChangeOfValues);	                \
    }										\
  }										\
  for(iPost = 0 ; iPost < List_Nbr(PE_L) ; iPost++)				\
     Destroy_PostElement(*(struct PostElement **)List_Pointer(PE_L, iPost));
  

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

  struct CutEdge       e[NBR_MAX_CUT];
  struct Element       Element ;
  struct PostElement * PE ;
  struct Value       * CumulativeValues ;
  List_T             * PE_L ;

  int     NbGeoElement, NbTimeStep, NbCut, * NumNodes ;
  int     iPost, iNode, iGeo, iCut, iEdge, iTime ;
  double  A, B, C, D, d1, d2, u, xcg, ycg, zcg ;
  double  x[3], y[3], z[3] ;

  GetDP_Begin("Pos_PrintOnSection");

  NbTimeStep = Pos_InitTimeSteps(PostSubOperation_P);

  PE_L = List_Create(10, 10, sizeof(struct PostElement *)) ;

  for(iCut = 0 ; iCut < NBR_MAX_CUT ; iCut++)
    e[iCut].Value = (struct Value*) Malloc(NbTimeStep*sizeof(struct Value)) ;    
    
  Format_PostHeader(PostSubOperation_P->Format, 
		    PostSubOperation_P->Iso, NbTimeStep, 
		    PostSubOperation_P->HarmonicToTime,
		    PostSubOperation_P->CombinationType, Order,
		    NCPQ_P?NCPQ_P->Name:NULL,
		    CPQ_P?CPQ_P->Name:NULL);

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

  switch(PostSubOperation_P->SubType) {

  case PRINT_ONSECTION_1D :
    Msg(GERROR, "Print on 1D cuts not done (use Print OnLine instead)");
    break;

  case PRINT_ONSECTION_2D :
    
    /* Ax+By+Cz+D=0  from  (x0,y0,z0),(x1,y1,z1),(x2,y2,z2) */
    
    x[0] = PostSubOperation_P->Case.OnSection.x[0] ;
    y[0] = PostSubOperation_P->Case.OnSection.y[0] ;
    z[0] = PostSubOperation_P->Case.OnSection.z[0] ;
    x[1] = PostSubOperation_P->Case.OnSection.x[1] ;
    y[1] = PostSubOperation_P->Case.OnSection.y[1] ;
    z[1] = PostSubOperation_P->Case.OnSection.z[1] ;
    x[2] = PostSubOperation_P->Case.OnSection.x[2] ;
    y[2] = PostSubOperation_P->Case.OnSection.y[2] ;
    z[2] = PostSubOperation_P->Case.OnSection.z[2] ;
    A =  (y[1]-y[0])*(z[2]-z[0]) - (z[1]-z[0])*(y[2]-y[0]) ;
    B = -(x[1]-x[0])*(z[2]-z[0]) + (z[1]-z[0])*(x[2]-x[0]) ;
    C =  (x[1]-x[0])*(y[2]-y[0]) - (y[1]-y[0])*(x[2]-x[0]) ;
    D = -A*x[0]-B*y[0]-C*z[0] ;

    /* Cut each element */
    
    NbGeoElement = Geo_GetNbrGeoElements() ;

    for(iGeo = 0 ; iGeo < NbGeoElement ; iGeo++) {
      if(InteractiveInterrupt) break;   
      Progress(iGeo, NbGeoElement, "Cut: ") ;
      Element.GeoElement = Geo_GetGeoElement(iGeo) ;
      Element.Num        = Element.GeoElement->Num ;
      Element.Type       = Element.GeoElement->Type ;
      Current.Region = Element.Region = Element.GeoElement->Region ;

      if((PostSubOperation_P->Dimension == _ALL && 
	  (Element.GeoElement->Type != POINT)) ||
	 (PostSubOperation_P->Dimension == _3D && 
	  (Element.GeoElement->Type & (TETRAHEDRON|HEXAHEDRON|PRISM|PYRAMID))) ||
	 (PostSubOperation_P->Dimension == _2D && 
	  (Element.GeoElement->Type & (TRIANGLE|QUADRANGLE))) ||
	 (PostSubOperation_P->Dimension == _1D && 
	  (Element.GeoElement->Type & LINE))){

	Get_NodesCoordinatesOfElement(&Element) ;

	if(Element.GeoElement->NbrEdges == 0)
	  Geo_CreateEdgesOfElement(Element.GeoElement) ;

	NbCut = 0;
	
	for(iEdge = 0 ; iEdge < Element.GeoElement->NbrEdges ; iEdge++){
	  NumNodes = Geo_GetNodesOfEdgeInElement(Element.GeoElement, iEdge) ;
	  e[NbCut].x[0] = Element.x[abs(NumNodes[0])-1] ;