pos_format.c

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

  double  TimeMH ;
  struct Value  TmpValue ;

  GetDP_Begin("Format_GmshParsed");

  if(Gmsh_StartNewView){
    Gmsh_StartNewView = 0 ;
    if(!TimeValue_L) TimeValue_L = List_Create(100,1000,sizeof(double));
    else List_Reset(TimeValue_L);
  }

  if (HarmonicToTime == 1)
    for(k = 0 ; k < NbHarmonic ; k++)
      List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
  else
    for(k = 0 ; k < HarmonicToTime ; k++)
      List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);	
  
  switch (Value[0].Type) {

  case SCALAR :

    if(TimeStep == 0){
      switch(Type){
      case POINT       : fprintf(PostStream, "SP("); break;
      case LINE        : fprintf(PostStream, "SL("); break;
      case TRIANGLE    : fprintf(PostStream, "ST("); break;
      case QUADRANGLE  : fprintf(PostStream, "SQ("); break;
      case TETRAHEDRON : fprintf(PostStream, "SS("); break;
      case HEXAHEDRON  : fprintf(PostStream, "SH("); break;
      case PRISM       : fprintf(PostStream, "SI("); break;
      case PYRAMID     : fprintf(PostStream, "SY("); break;
      }    
      for(i = 0 ; i < NbrNodes ; i++){
	if(i) fprintf(PostStream, ",");
	fprintf(PostStream, "%.8g,%.8g,%.8g", x[i], y[i], z[i]);
      }
      fprintf(PostStream, "){");
    }

    if (HarmonicToTime == 1) {
      for(k = 0 ; k < NbHarmonic ; k++) {
	if(k || TimeStep) fprintf(PostStream, ",");
	for(i = 0 ; i < NbrNodes ; i++){
	  if(i) fprintf(PostStream, ",");
	  fprintf(PostStream, "%.8g", Value[i].Val[MAX_DIM*k]);
	}
      }
    }
    else {
      for(k = 0 ; k < HarmonicToTime ; k++){
	if(k || TimeStep) fprintf(PostStream, ",");
	for(i = 0 ; i < NbrNodes ; i++){
	  if(i) fprintf(PostStream, ",");
	  F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
	  fprintf(PostStream, "%.8g", TmpValue.Val[0]);
	}
      }
    }

    if(TimeStep == NbTimeStep-1){
      fprintf(PostStream, "};\n") ;
    }
    break ;

  case VECTOR :

    if(TimeStep == 0){
      switch(Type){
      case POINT       : fprintf(PostStream, "VP("); break;
      case LINE        : fprintf(PostStream, "VL("); break;
      case TRIANGLE    : fprintf(PostStream, "VT("); break;
      case QUADRANGLE  : fprintf(PostStream, "VQ("); break;
      case TETRAHEDRON : fprintf(PostStream, "VS("); break;
      case HEXAHEDRON  : fprintf(PostStream, "VH("); break;
      case PRISM       : fprintf(PostStream, "VI("); break;
      case PYRAMID     : fprintf(PostStream, "VY("); break;
      }
      for(i = 0 ; i < NbrNodes ; i++){
	if(i) fprintf(PostStream, ",");
	fprintf(PostStream, "%.8g,%.8g,%.8g", x[i], y[i], z[i]);
      }
      fprintf(PostStream, "){");
    }
    
    if (HarmonicToTime == 1) {
      for(k = 0 ; k < NbHarmonic ; k++) {
	if(k || TimeStep) fprintf(PostStream, ",");
	for(i = 0 ; i < NbrNodes ; i++){
	  if(i) fprintf(PostStream, ",");
	  for(j = 0 ; j < 3 ; j++){
	    if(j) fprintf(PostStream, ",");
	    fprintf(PostStream, "%.8g", Value[i].Val[MAX_DIM*k+j]);
	  }
	}
      }
    }
    else {
      for(k = 0 ; k < HarmonicToTime ; k++){
	if(k || TimeStep) fprintf(PostStream, ",");
	for(i = 0 ; i < NbrNodes ; i++){
	  if(i) fprintf(PostStream, ",");
	  F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
	  for(j = 0 ; j < 3 ; j++){
	    if(j) fprintf(PostStream, ",");
	    fprintf(PostStream, "%.8g", TmpValue.Val[j]);
	  }
	}
      }
    }

    if(TimeStep == NbTimeStep-1){
      fprintf(PostStream, "};\n") ;
    }
    break ;

  case TENSOR_DIAG :
  case TENSOR_SYM :
  case TENSOR :
    Msg(GERROR, "Gmsh parsed format not done for Tensors");
    break;
    
  }

  GetDP_End ;
}


void  Format_Gmsh(double Time, int TimeStep, int NbTimeStep, int NbHarmonic, 
		  int HarmonicToTime, int Type, int NbrNodes, 
		  double *x, double *y, double *z, struct Value *Value){
  int            i,j,k;
  double         TimeMH ;
  struct Value   TmpValue ;
  static List_T *Current_L ;

  GetDP_Begin("Format_Gmsh");

  if(Gmsh_StartNewView){
    Gmsh_StartNewView = 0 ;
    NbSP = NbVP = NbTP = NbSL = NbVL = NbTL = 0;
    NbST = NbVT = NbTT = NbSQ = NbVQ = NbTQ = 0;
    NbSS = NbVS = NbTS = NbSH = NbVH = NbTH = 0;
    NbSI = NbVI = NbTI = NbSY = NbVY = NbTY = 0;
    NbT2 = 0;
    if(!SP) SP = List_Create(100,1000,sizeof(double)); else List_Reset(SP);
    if(!VP) VP = List_Create(100,1000,sizeof(double)); else List_Reset(VP);
    if(!TP) TP = List_Create(100,1000,sizeof(double)); else List_Reset(TP);
    if(!SL) SL = List_Create(100,1000,sizeof(double)); else List_Reset(SL);
    if(!VL) VL = List_Create(100,1000,sizeof(double)); else List_Reset(VL);
    if(!TL) TL = List_Create(100,1000,sizeof(double)); else List_Reset(TL);
    if(!ST) ST = List_Create(100,1000,sizeof(double)); else List_Reset(ST);
    if(!VT) VT = List_Create(100,1000,sizeof(double)); else List_Reset(VT);
    if(!TT) TT = List_Create(100,1000,sizeof(double)); else List_Reset(TT);
    if(!SQ) SQ = List_Create(100,1000,sizeof(double)); else List_Reset(SQ);
    if(!VQ) VQ = List_Create(100,1000,sizeof(double)); else List_Reset(VQ);
    if(!TQ) TQ = List_Create(100,1000,sizeof(double)); else List_Reset(TQ);
    if(!SS) SS = List_Create(100,1000,sizeof(double)); else List_Reset(SS);
    if(!VS) VS = List_Create(100,1000,sizeof(double)); else List_Reset(VS);
    if(!TS) TS = List_Create(100,1000,sizeof(double)); else List_Reset(TS);
    if(!SH) SH = List_Create(100,1000,sizeof(double)); else List_Reset(SH);
    if(!VH) VH = List_Create(100,1000,sizeof(double)); else List_Reset(VH);
    if(!TH) TH = List_Create(100,1000,sizeof(double)); else List_Reset(TH);
    if(!SI) SI = List_Create(100,1000,sizeof(double)); else List_Reset(SI);
    if(!VI) VI = List_Create(100,1000,sizeof(double)); else List_Reset(VI);
    if(!TI) TI = List_Create(100,1000,sizeof(double)); else List_Reset(TI);
    if(!SY) SY = List_Create(100,1000,sizeof(double)); else List_Reset(SY);
    if(!VY) VY = List_Create(100,1000,sizeof(double)); else List_Reset(VY);
    if(!TY) TY = List_Create(100,1000,sizeof(double)); else List_Reset(TY);
    if(!T2D) T2D = List_Create(100,1000,sizeof(double)); else List_Reset(T2D);
    if(!T2C) T2C = List_Create(100,1000,sizeof(char)); else List_Reset(T2C);
    if(!TimeValue_L) TimeValue_L = List_Create(100,1000,sizeof(double));
    else List_Reset(TimeValue_L);
  }

  switch (Value[0].Type) {

  case SCALAR :

    if(TimeStep == 0){
      switch(Type){
      case POINT       : Current_L = SP ; NbSP++ ; break ;
      case LINE        : Current_L = SL ; NbSL++ ; break ;
      case TRIANGLE    : Current_L = ST ; NbST++ ; break ;
      case QUADRANGLE  : Current_L = SQ ; NbSQ++ ; break ;
      case TETRAHEDRON : Current_L = SS ; NbSS++ ; break ;
      case HEXAHEDRON  : Current_L = SH ; NbSH++ ; break ;
      case PRISM       : Current_L = SI ; NbSI++ ; break ;
      case PYRAMID     : Current_L = SY ; NbSY++ ; break ;
      }
      for(i = 0 ; i < NbrNodes ; i++) List_Add(Current_L, &x[i]);
      for(i = 0 ; i < NbrNodes ; i++) List_Add(Current_L, &y[i]);
      for(i = 0 ; i < NbrNodes ; i++) List_Add(Current_L, &z[i]);
    }
    if (HarmonicToTime == 1)
      for(k = 0 ; k < NbHarmonic ; k++){
	List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
	for(i = 0 ; i < NbrNodes ; i++)
	  List_Add(Current_L, &Value[i].Val[MAX_DIM*k]);
      }
    else
      for(k = 0 ; k < HarmonicToTime ; k++){
	List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);	
	for(i = 0 ; i < NbrNodes ; i++){
	  F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
	  List_Add(Current_L, &TmpValue.Val[0]);
	}
      }
    break ;

  case VECTOR :

    if(TimeStep == 0){
      switch(Type){
      case POINT       : Current_L = VP ; NbVP++ ; break ;
      case LINE        : Current_L = VL ; NbVL++ ; break ;
      case TRIANGLE    : Current_L = VT ; NbVT++ ; break ;
      case QUADRANGLE  : Current_L = VQ ; NbVQ++ ; break ;
      case TETRAHEDRON : Current_L = VS ; NbVS++ ; break ;
      case HEXAHEDRON  : Current_L = VH ; NbVH++ ; break ;
      case PRISM       : Current_L = VI ; NbVI++ ; break ;
      case PYRAMID     : Current_L = VY ; NbVY++ ; break ;
      }
      for(i = 0 ; i < NbrNodes ; i++) List_Add(Current_L, &x[i]);
      for(i = 0 ; i < NbrNodes ; i++) List_Add(Current_L, &y[i]);
      for(i = 0 ; i < NbrNodes ; i++) List_Add(Current_L, &z[i]);
    }
    if (HarmonicToTime == 1)
      for(k = 0 ; k < NbHarmonic ; k++){
	List_Put(TimeValue_L, NbHarmonic*TimeStep+k, &Time);
	for(i = 0 ; i < NbrNodes ; i++)
	  for(j = 0 ; j < 3 ; j++)
	    List_Add(Current_L, &Value[i].Val[MAX_DIM*k+j]);
      }
    else
      for(k = 0 ; k < HarmonicToTime ; k++){
	List_Put(TimeValue_L, HarmonicToTime*TimeStep+k, &Time);
	for(i = 0 ; i < NbrNodes ; i++){
	  F_MHToTime0(k+i, &Value[i], &TmpValue, k, HarmonicToTime, &TimeMH) ;
	  for(j = 0 ; j < 3 ; j++)
	    List_Add(Current_L, &TmpValue.Val[j]);
	}
      }
    break ;

  case TENSOR_DIAG :
  case TENSOR_SYM :
  case TENSOR :
    Msg(GERROR, "Gmsh format not done for Tensors");
    break;
    
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  F o r m a t _ G n u p l o t                                             */
/* ------------------------------------------------------------------------ */

int Get_GmshElementType(int Type){

  GetDP_Begin("Get_GmshElementType");

  switch(Type){
  case POINT :         GetDP_Return(15) ; 
  case LINE :	       GetDP_Return(1)  ; 
  case TRIANGLE :      GetDP_Return(2)  ; 
  case QUADRANGLE :    GetDP_Return(3)  ; 
  case TETRAHEDRON :   GetDP_Return(4)  ;    
  case HEXAHEDRON :    GetDP_Return(5)  ; 
  case PRISM :	       GetDP_Return(6)  ; 
  case PYRAMID :       GetDP_Return(7)  ; 
  case LINE_2 :	       GetDP_Return(8)  ; 
  case TRIANGLE_2 :    GetDP_Return(9)  ; 
  case QUADRANGLE_2 :  GetDP_Return(10) ; 
  case TETRAHEDRON_2 : GetDP_Return(11) ;      
  case HEXAHEDRON_2 :  GetDP_Return(12) ; 
  case PRISM_2 :       GetDP_Return(13) ; 
  case PYRAMID_2 :     GetDP_Return(14) ; 
  default : 
    Msg(GERROR, "Unknown type of element in Gmsh format") ;
    GetDP_Return(-1) ;
  }
}

void  Format_Gnuplot(int Format, double Time, int TimeStep, int NbrTimeSteps,
		     int NbrHarmonics, int HarmonicToTime, 
		     int ElementType, int NumElement, int NbrNodes,
		     double *x, double *y, double *z, double *Dummy,
		     struct Value *Value){
  static int  Size, TmpIndex ;
  static double  * TmpValues ;
  int         i, j, k, t, i2, k2 ;
  double      TimeMH ;
  struct Value  TmpValue ;

  GetDP_Begin("Format_Gnuplot");

  if(TimeStep == 0){
    switch(Value->Type){
    case SCALAR      : Size = 1 ; break ;
    case VECTOR      : Size = 3 ; break ;
    case TENSOR_DIAG : Size = 3 ; break ;
    case TENSOR_SYM  : Size = 6 ; break ;
    case TENSOR      : Size = 9 ; break ;
    }
    TmpValues = (double*) Malloc(NbrTimeSteps*NbrNodes*NbrHarmonics*Size*sizeof(double));
    TmpIndex = 0;
  }

  for(i = 0 ; i < NbrNodes ; i++)
    for(k = 0 ; k < NbrHarmonics ; k++)
      for(j = 0 ; j < Size ; j++)
	TmpValues[TmpIndex++] = Value[i].Val[MAX_DIM*k+j];

  if(TimeStep == NbrTimeSteps-1){

    for(i = 0 ; i <= NbrNodes ; i++){ /* New line for each node, closed loop for tri/qua */
      
      if(i != NbrNodes)
	i2 = i ;
      else{
	if(NbrNodes < 3) break ;
	else i2 = 0 ;
      }

      fprintf(PostStream, "%d %d ", Get_GmshElementType(ElementType), NumElement);
      fprintf(PostStream, " %.16g %.16g %.16g ", x[i2], y[i2], z[i2]);
      if(Dummy){
	if(Dummy[3]<0){
	  if(!i)
	    fprintf(PostStream, " %.16g %.16g 0 ", Dummy[0], Dummy[2]);
	  else
	    fprintf(PostStream, " %.16g %.16g 0 ", Dummy[1], Dummy[2]);
	}
	else
	  fprintf(PostStream, " %.16g %.16g %.16g ", Dummy[0], Dummy[1], Dummy[2]);
      }
      else
	fprintf(PostStream, " 0 0 0 ");

      for(t = 0 ; t < NbrTimeSteps ; t++){

	if (HarmonicToTime == 1) {
	  for(k = 0 ; k < NbrHarmonics ; k++) {
	    for(j = 0 ; j < Size ; j++){
	      fprintf(PostStream, " %.16g", 
		      TmpValues[ t*NbrNodes*NbrHarmonics*Size
			       + i2*NbrHarmonics*Size
			       + k*Size
			       + j ]);
	    }
	    fprintf(PostStream, " ");
	  }
	}
	else {
	  TmpValue.Type = Value->Type ;
	  for(k = 0 ; k < HarmonicToTime ; k++){

	    for(k2 = 0 ; k2 < NbrHarmonics ; k2++)
	      for(j = 0 ; j < Size ; j++)
		TmpValue.Val[MAX_DIM*k2+j] =
		  TmpValues[ t*NbrNodes*NbrHarmonics*Size
			   + i2*NbrHarmonics*Size
			   + k2*Size
			   + j ] ;

	    F_MHToTime0(k, &TmpValue, &TmpValue, k, HarmonicToTime, &TimeMH) ;
	    for(j = 0 ; j < Size ; j++)
	      fprintf(PostStream, "%.16g", TmpValue.Val[0]);
	    fprintf(PostStream, " ");
	  }
	}
	fprintf(PostStream, " ");

      } /* for t */
      fprintf(PostStream, "\n");      

    } /* for i */