f_type.c

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

      V->Val[MAX_DIM*k+7] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+7] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+8] ; Im = A->Val[MAX_DIM*(k+1)+8] ;
      V->Val[MAX_DIM*k+8] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+8] = atan2(Im,Re);   
    }
    break;

  default :
    Msg(GERROR, "Unknown type of arguments in function 'Cart2Pol'");
    break;
  }
  
  V->Type = A->Type ;

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  Create 1 Vector from 3 Scalar                                           */
/* ------------------------------------------------------------------------ */

void  F_Vector (F_ARG) {
  int k ;

  GetDP_Begin("F_Vector");

  if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
    Msg(GERROR, "Non scalar argument(s) for function 'Vector'");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = (A  )->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
  }
  V->Type = VECTOR ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Create 1 Tensor from 9 Scalar                                           */
/* ------------------------------------------------------------------------ */

void  F_Tensor (F_ARG) {
  int k ;

  GetDP_Begin("F_Tensor");

  if(  (A)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR ||
     (A+3)->Type != SCALAR || (A+4)->Type != SCALAR || (A+5)->Type != SCALAR ||
     (A+6)->Type != SCALAR || (A+7)->Type != SCALAR || (A+8)->Type != SCALAR )
    Msg(GERROR, "Non scalar argument(s) for function 'Tensor'");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = (A  )->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+3] = (A+3)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+4] = (A+4)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+5] = (A+5)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+6] = (A+6)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+7] = (A+7)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+8] = (A+8)->Val[MAX_DIM*k] ;
  }
  V->Type = TENSOR ;

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  Create 1 Symmetric Tensor from 6 Scalar                                 */
/* ------------------------------------------------------------------------ */

void  F_TensorSym (F_ARG) {
  int k ;

  GetDP_Begin("F_TensorSym");

  if(  (A)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR ||
     (A+3)->Type != SCALAR || (A+4)->Type != SCALAR || (A+5)->Type != SCALAR )
    Msg(GERROR, "Non scalar argument(s) for function 'TensorSym'");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = (A  )->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+3] = (A+3)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+4] = (A+4)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+5] = (A+5)->Val[MAX_DIM*k] ;
  }
  V->Type = TENSOR_SYM ;

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  Create 1 Diagonal Tensor from 3 Scalar                                  */
/* ------------------------------------------------------------------------ */

void  F_TensorDiag (F_ARG) {
  int k ;

  GetDP_Begin("F_TensorDiag");

  if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
    Msg(GERROR, "Non scalar argument(s) for function 'TensorDiag'");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] =     A->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+1] = (A+1)->Val[MAX_DIM*k] ;
    V->Val[MAX_DIM*k+2] = (A+2)->Val[MAX_DIM*k] ;
  }
  V->Type = TENSOR_DIAG ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Create 1 Tensor from 3 Vector                                           */
/* ------------------------------------------------------------------------ */

void  F_TensorV (F_ARG) {
  int k ;

  GetDP_Begin("F_TensorV");

  if((A)->Type != VECTOR || (A+1)->Type != VECTOR || (A+2)->Type != VECTOR)
    Msg(GERROR, "Non scalar argument(s) for function 'TensorV'");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = (A  )->Val[MAX_DIM*k  ] ;
    V->Val[MAX_DIM*k+1] = (A  )->Val[MAX_DIM*k+1] ;
    V->Val[MAX_DIM*k+2] = (A  )->Val[MAX_DIM*k+2] ;
    V->Val[MAX_DIM*k+3] = (A+1)->Val[MAX_DIM*k  ] ;
    V->Val[MAX_DIM*k+4] = (A+1)->Val[MAX_DIM*k+1] ;
    V->Val[MAX_DIM*k+5] = (A+1)->Val[MAX_DIM*k+2] ;
    V->Val[MAX_DIM*k+6] = (A+2)->Val[MAX_DIM*k  ] ;
    V->Val[MAX_DIM*k+7] = (A+2)->Val[MAX_DIM*k+1] ;
    V->Val[MAX_DIM*k+8] = (A+2)->Val[MAX_DIM*k+2] ;
  }
  V->Type = TENSOR ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Dyadic product                                                          */
/* ------------------------------------------------------------------------ */

void  F_SquDyadicProduct (F_ARG) {
  int k ;
  double t11, t12, t13, t22, t23, t33 ;

  GetDP_Begin("F_SquDyadicProduct");

  if (A->Type != VECTOR)
    Msg(GERROR, "Non vector argument for function 'TensorDyadic'");

  t11 = SQU(A->Val[0]) ;  t22 = SQU(A->Val[1]) ;  t33 = SQU(A->Val[2]) ;
  t12 = A->Val[0] * A->Val[1] ;  t13 = A->Val[0] * A->Val[2] ;
  t23 = A->Val[1] * A->Val[2] ;

  V->Val[0] = t11 ;  V->Val[1] = t12 ;  V->Val[2] = t13 ;
  V->Val[3] = t22 ;  V->Val[4] = t23 ;  V->Val[5] = t33 ;

  /* Attention : a revoir */
  if (Current.NbrHar > 1) {
    V->Val[MAX_DIM  ] = V->Val[MAX_DIM+1] = V->Val[MAX_DIM+2] =
      V->Val[MAX_DIM+3] = V->Val[MAX_DIM+4] = V->Val[MAX_DIM+5] = 0. ;
    for (k = 2 ; k < Current.NbrHar ; k++) {
      V->Val[MAX_DIM*k  ] = V->Val[MAX_DIM*k+1] = V->Val[MAX_DIM*k+2] =
	V->Val[MAX_DIM*k+3] = V->Val[MAX_DIM*k+4] = V->Val[MAX_DIM*k+5] = 0. ;
    }
  }

  V->Type = TENSOR_SYM ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Get Vector Components                                                   */
/* ------------------------------------------------------------------------ */

#define get_comp_vector(index,string)					\
  int k ;								\
									\
  GetDP_Begin("F_" string);						\
									\
  if(A->Type != VECTOR)							\
    Msg(GERROR, "Non vector argument for function '" string "'");	\
									\
  for (k = 0 ; k < Current.NbrHar ; k++) {				\
    V->Val[MAX_DIM*k  ] = A->Val[MAX_DIM*k+index] ;			\
  }									\
  V->Type = SCALAR ;							\
									\
  GetDP_End ;
    
void  F_CompX (F_ARG) { get_comp_vector(0,"CompX") }
void  F_CompY (F_ARG) { get_comp_vector(1,"CompY") }
void  F_CompZ (F_ARG) { get_comp_vector(2,"CompZ") }

#undef get_comp_vector


/* ------------------------------------------------------------------------ */
/*  Get Tensor Components                                                   */
/* ------------------------------------------------------------------------ */

#define get_comp_tensor(i,is,id,string)							\
  int k ;										\
											\
  GetDP_Begin("F_" string);								\
											\
  switch(A->Type) {									\
  case TENSOR :										\
    for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k+(i)] ;	\
    break ;										\
  case TENSOR_SYM :									\
    for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k+(is)] ;	\
    break ;										\
  case TENSOR_DIAG :									\
    if(id >= 0)										\
      for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = A->Val[MAX_DIM*k+(id)] ;	\
    else										\
      for (k=0; k<Current.NbrHar; k++) V->Val[MAX_DIM*k] = 0.;				\
    break ;										\
  default :										\
    Msg(GERROR, "Non tensor argument for function '" string "'");			\
    break;										\
  }											\
  V->Type = SCALAR ;									\
  											\
  GetDP_End ;
  
void  F_CompXX (F_ARG) { get_comp_tensor(0,0, 0,"CompXX") }
void  F_CompXY (F_ARG) { get_comp_tensor(1,1,-1,"CompXY") }
void  F_CompXZ (F_ARG) { get_comp_tensor(2,2,-1,"CompXZ") }
void  F_CompYX (F_ARG) { get_comp_tensor(3,1,-1,"CompYX") }
void  F_CompYY (F_ARG) { get_comp_tensor(4,3, 1,"CompYY") }
void  F_CompYZ (F_ARG) { get_comp_tensor(5,4,-1,"CompYZ") }
void  F_CompZX (F_ARG) { get_comp_tensor(6,2,-1,"CompZX") }
void  F_CompZY (F_ARG) { get_comp_tensor(7,4,-1,"CompZY") }
void  F_CompZZ (F_ARG) { get_comp_tensor(8,5, 2,"CompZZ") }

#undef get_comp_tensor




/* ------------------------------------------------------------------------ */
/*  Get Tensor for transformation of vector                                 */ 
/*  from cartesian to spherical coordinate system                           */
/* ------------------------------------------------------------------------ */

void  F_Cart2Sph (F_ARG) {
  int k ;
  double theta, phi ;

  GetDP_Begin("F_Cart2Sph");
  
  if((A)->Type != VECTOR)
     Msg(GERROR, "Vector argument required for Function 'Cart2Sph'");
     
  /* Warning! This is the physic's convention. For the math
     convention, switch theta and phi. */

  theta = atan2( sqrt(SQU(A->Val[0])+ SQU(A->Val[1])) ,  A->Val[2] ) ;   
  phi = atan2( A->Val[1] , A->Val[0] ) ; 
    
  /* r basis vector */
  V->Val[0] = sin(theta) * cos(phi) ; 
  V->Val[1] = sin(theta) * sin(phi) ; 
  V->Val[2] = cos(theta) ;

  /* theta basis vector */
  V->Val[3] = cos(theta) * cos(phi) ; 
  V->Val[4] = cos(theta) * sin(phi) ; 
  V->Val[5] = - sin(theta) ; 

  /* phi basis vector */
  V->Val[6] = - sin(phi) ; 
  V->Val[7] = cos(phi) ; 
  V->Val[8] = 0. ;
    
  for (k = 0 ; k < Current.NbrHar ; k+=2) {
    V->Val[MAX_DIM*k  ] = V->Val[0] ;
    V->Val[MAX_DIM*k+1] = V->Val[1] ;
    V->Val[MAX_DIM*k+2] = V->Val[2] ;
    V->Val[MAX_DIM*k+3] = V->Val[3] ;
    V->Val[MAX_DIM*k+4] = V->Val[4] ;
    V->Val[MAX_DIM*k+5] = V->Val[5] ;
    V->Val[MAX_DIM*k+6] = V->Val[6] ;
    V->Val[MAX_DIM*k+7] = V->Val[7] ;
    V->Val[MAX_DIM*k+8] = V->Val[8] ;
    V->Val[MAX_DIM*(k+1)  ] = 0. ;
    V->Val[MAX_DIM*(k+1)+1] = 0. ;
    V->Val[MAX_DIM*(k+1)+2] = 0. ;
    V->Val[MAX_DIM*(k+1)+3] = 0. ;
    V->Val[MAX_DIM*(k+1)+4] = 0. ;
    V->Val[MAX_DIM*(k+1)+5] = 0. ;
    V->Val[MAX_DIM*(k+1)+6] = 0. ;
    V->Val[MAX_DIM*(k+1)+7] = 0. ;
    V->Val[MAX_DIM*(k+1)+8] = 0. ;
  }
  V->Type = TENSOR ;

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  Get Tensor for transformation of vector                                 */ 
/*  from cartesian to cylindric coordinate system                           */
/*  vector              ->  Cart2Cyl[XYZ[]] * vector                        */
/*  (x,y,z)-components  ->  (radial, tangential, axial)-components          */
/* ------------------------------------------------------------------------ */

void  F_Cart2Cyl (F_ARG) {
  int k ;
  double theta ;

  GetDP_Begin("F_Cart2Cyl");
  
  if((A)->Type != VECTOR)
     Msg(GERROR, "Vector argument required for Function 'Cart2Cyl'");
     
  theta = atan2(A->Val[1] , A->Val[0]) ;   
    
  V->Val[0] =  cos(theta) ; 
  V->Val[1] =  sin(theta) ; 
  V->Val[2] =  0 ;
  V->Val[3] = -sin(theta) ; 
  V->Val[4] =  cos(theta) ; 
  V->Val[5] =  0 ; 
  V->Val[6] =  0 ; 
  V->Val[7] =  0 ; 
  V->Val[8] =  1. ;
    
  for (k = 0 ; k < Current.NbrHar ; k+=2) {
    V->Val[MAX_DIM*k  ] = V->Val[0] ;
    V->Val[MAX_DIM*k+1] = V->Val[1] ;
    V->Val[MAX_DIM*k+2] = V->Val[2] ;
    V->Val[MAX_DIM*k+3] = V->Val[3] ;
    V->Val[MAX_DIM*k+4] = V->Val[4] ;
    V->Val[MAX_DIM*k+5] = V->Val[5] ;
    V->Val[MAX_DIM*k+6] = V->Val[6] ;
    V->Val[MAX_DIM*k+7] = V->Val[7] ;
    V->Val[MAX_DIM*k+8] = V->Val[8] ;
    V->Val[MAX_DIM*(k+1)  ] = 0. ;
    V->Val[MAX_DIM*(k+1)+1] = 0. ;
    V->Val[MAX_DIM*(k+1)+2] = 0. ;
    V->Val[MAX_DIM*(k+1)+3] = 0. ;
    V->Val[MAX_DIM*(k+1)+4] = 0. ;
    V->Val[MAX_DIM*(k+1)+5] = 0. ;
    V->Val[MAX_DIM*(k+1)+6] = 0. ;
    V->Val[MAX_DIM*(k+1)+7] = 0. ;
    V->Val[MAX_DIM*(k+1)+8] = 0. ;
  }
  V->Type = TENSOR ;

  GetDP_End ;
}




/* ------------------------------------------------------------------------ */
/*  U n i t V e c t o r X, Y, Z                                             */
/* ------------------------------------------------------------------------ */

void  F_UnitVectorX (F_ARG) {
  int k ;

  GetDP_Begin("F_UnitVectorX");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = (k)? 0.:1. ;
    V->Val[MAX_DIM*k+1] = 0. ;
    V->Val[MAX_DIM*k+2] = 0. ;
  }
  V->Type = VECTOR ;

  GetDP_End ;
}

void  F_UnitVectorY (F_ARG) {
  int k ;

  GetDP_Begin("F_UnitVectorY");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = 0. ;
    V->Val[MAX_DIM*k+1] = (k)? 0.:1. ;
    V->Val[MAX_DIM*k+2] = 0. ;
  }
  V->Type = VECTOR ;

  GetDP_End ;
}

void  F_UnitVectorZ (F_ARG) {
  int k ;

  GetDP_Begin("F_UnitVectorZ");

  for (k = 0 ; k < Current.NbrHar ; k++) {
    V->Val[MAX_DIM*k  ] = 0. ;
    V->Val[MAX_DIM*k+1] = 0. ;
    V->Val[MAX_DIM*k+2] = (k)? 0.:1. ;
  }
  V->Type = VECTOR ;

  GetDP_End ;
}

#undef F_ARG