cal_value.c

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

		   V1->Val[MAX_DIM*k+1] != V2->Val[MAX_DIM*k+1] ||
		   V1->Val[MAX_DIM*k+2] != V2->Val[MAX_DIM*k+2]) ;
    }
    R->Type = SCALAR ;
  }
  else if ( (V1->Type == TENSOR_SYM) && (V2->Type == TENSOR_SYM) ) {
    R->Val[0] = (V1->Val[0] != V2->Val[0] ||
		 V1->Val[1] != V2->Val[1] ||
		 V1->Val[2] != V2->Val[2] ||
		 V1->Val[3] != V2->Val[3] ||
		 V1->Val[4] != V2->Val[4] ||
		 V1->Val[5] != V2->Val[5]) ;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      if(R->Val[0]) break;
      R->Val[0] = (V1->Val[MAX_DIM*k  ] != V2->Val[MAX_DIM*k  ] ||
		   V1->Val[MAX_DIM*k+1] != V2->Val[MAX_DIM*k+1] ||
		   V1->Val[MAX_DIM*k+2] != V2->Val[MAX_DIM*k+2] ||
		   V1->Val[MAX_DIM*k+3] != V2->Val[MAX_DIM*k+3] ||
		   V1->Val[MAX_DIM*k+4] != V2->Val[MAX_DIM*k+4] ||
		   V1->Val[MAX_DIM*k+5] != V2->Val[MAX_DIM*k+5]) ;
    }
    R->Type = SCALAR ;
  }
  else if ( (V1->Type == TENSOR) && (V2->Type == TENSOR) ) {
    R->Val[0] = (V1->Val[0] != V2->Val[0] ||
		 V1->Val[1] != V2->Val[1] ||
		 V1->Val[2] != V2->Val[2] ||
		 V1->Val[3] != V2->Val[3] ||
		 V1->Val[4] != V2->Val[4] ||
		 V1->Val[5] != V2->Val[5] ||
		 V1->Val[6] != V2->Val[6] ||
		 V1->Val[7] != V2->Val[7] ||
		 V1->Val[8] != V2->Val[8]) ;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      if(R->Val[0]) break;
      R->Val[0] = (V1->Val[MAX_DIM*k  ] != V2->Val[MAX_DIM*k  ] ||
		   V1->Val[MAX_DIM*k+1] != V2->Val[MAX_DIM*k+1] ||
		   V1->Val[MAX_DIM*k+2] != V2->Val[MAX_DIM*k+2] ||
		   V1->Val[MAX_DIM*k+3] != V2->Val[MAX_DIM*k+3] ||
		   V1->Val[MAX_DIM*k+4] != V2->Val[MAX_DIM*k+4] ||
		   V1->Val[MAX_DIM*k+5] != V2->Val[MAX_DIM*k+5] ||
		   V1->Val[MAX_DIM*k+6] != V2->Val[MAX_DIM*k+6] ||
		   V1->Val[MAX_DIM*k+7] != V2->Val[MAX_DIM*k+7] ||
		   V1->Val[MAX_DIM*k+8] != V2->Val[MAX_DIM*k+8]) ;
    }
    R->Type = SCALAR ;
  }
  else {
    Msg(GERROR, "Comparison of different quantities: %s != %s",
	Get_StringForDefine(Field_Type, V1->Type),
	Get_StringForDefine(Field_Type, V2->Type));
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ 
   R <- V1 ~= V2 
   ------------------------------------------------------------------------ */

void  Cal_ApproxEqualValue (struct Value * V1, struct Value * V2, struct Value * R) {

  GetDP_Begin("Cal_ApproxEqualValue");

  if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
    R->Val[0] = (fabs(V1->Val[0] - V2->Val[0]) < 1.e-10) ;
    R->Type = SCALAR ;
  }
  else {
    Msg(GERROR, "Comparison of non scalar quantities: %s ~= %s",
	Get_StringForDefine(Field_Type, V1->Type),
	Get_StringForDefine(Field_Type, V2->Type));
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ 
   R <- V1 && V2 
   ------------------------------------------------------------------------ */

void  Cal_AndValue (struct Value * V1, struct Value * V2, struct Value * R) {

  GetDP_Begin("Cal_AndValue");

  if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
    R->Val[0] = (V1->Val[0] && V2->Val[0]) ;
    R->Type = SCALAR ;
  }
  else {
    Msg(GERROR, "And of non scalar quantities: %s && %s",
	Get_StringForDefine(Field_Type, V1->Type),
	Get_StringForDefine(Field_Type, V2->Type));
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- V1 || V2 
   ------------------------------------------------------------------------ */

void  Cal_OrValue (struct Value * V1, struct Value * V2, struct Value * R) {

  GetDP_Begin("Cal_OrValue");

  if ( (V1->Type == SCALAR) && (V2->Type == SCALAR) ) {
    R->Val[0] = (V1->Val[0] || V2->Val[0]) ;
    R->Type = SCALAR ;
  }
  else {
    Msg(GERROR, "Or of non scalar quantities: %s || %s",
	Get_StringForDefine(Field_Type, V1->Type),
	Get_StringForDefine(Field_Type, V2->Type));
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- -R 
   ------------------------------------------------------------------------ */

void  Cal_NegValue (struct Value * R) {
  int k ;

  GetDP_Begin("Cal_NegValue");

  switch(R->Type) {
  case SCALAR :
    for (k = 0 ; k < Current.NbrHar ; k++){
      R->Val[MAX_DIM*k] = -R->Val[MAX_DIM*k] ;
    }
    break;
  case VECTOR :
  case TENSOR_DIAG :
    for (k = 0 ; k < Current.NbrHar ; k++){
      R->Val[MAX_DIM*k  ] = -R->Val[MAX_DIM*k  ] ;
      R->Val[MAX_DIM*k+1] = -R->Val[MAX_DIM*k+1] ;
      R->Val[MAX_DIM*k+2] = -R->Val[MAX_DIM*k+2] ;
    }
    break;
  case TENSOR_SYM :
    for (k = 0 ; k < Current.NbrHar ; k++){
      R->Val[MAX_DIM*k  ] = -R->Val[MAX_DIM*k  ] ;
      R->Val[MAX_DIM*k+1] = -R->Val[MAX_DIM*k+1] ;
      R->Val[MAX_DIM*k+2] = -R->Val[MAX_DIM*k+2] ;
      R->Val[MAX_DIM*k+3] = -R->Val[MAX_DIM*k+3] ;
      R->Val[MAX_DIM*k+4] = -R->Val[MAX_DIM*k+4] ;
      R->Val[MAX_DIM*k+5] = -R->Val[MAX_DIM*k+5] ;
    }
    break;
  case TENSOR :
    for (k = 0 ; k < Current.NbrHar ; k++){
      R->Val[MAX_DIM*k  ] = -R->Val[MAX_DIM*k  ] ;
      R->Val[MAX_DIM*k+1] = -R->Val[MAX_DIM*k+1] ;
      R->Val[MAX_DIM*k+2] = -R->Val[MAX_DIM*k+2] ;
      R->Val[MAX_DIM*k+3] = -R->Val[MAX_DIM*k+3] ;
      R->Val[MAX_DIM*k+4] = -R->Val[MAX_DIM*k+4] ;
      R->Val[MAX_DIM*k+5] = -R->Val[MAX_DIM*k+5] ;
      R->Val[MAX_DIM*k+6] = -R->Val[MAX_DIM*k+6] ;
      R->Val[MAX_DIM*k+7] = -R->Val[MAX_DIM*k+7] ;
      R->Val[MAX_DIM*k+8] = -R->Val[MAX_DIM*k+8] ;
    }
    break;
  default :
    Msg(GERROR, "Wrong argument type for Operator (-)");
    break;
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- !R 
   ------------------------------------------------------------------------ */

void  Cal_NotValue (struct Value * R) {

  GetDP_Begin("Cal_NotValue");

  if (R->Type == SCALAR){
    R->Val[0] = !R->Val[0] ;
  }
  else {
    Msg(GERROR, "Negation of non scalar quantity: ! %s",
	Get_StringForDefine(Field_Type, R->Type));
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- V1^T 
   ------------------------------------------------------------------------ */

void Cal_TransposeValue(struct Value *V1, struct Value *R){
  int     k;

  GetDP_Begin("Cal_TransposeValue");

  switch(V1->Type){

  case TENSOR_DIAG :
  case TENSOR_SYM :
    Cal_CopyValue(V1,R);
    break;

  case TENSOR :
    if(Current.NbrHar==1){      
      R->Val[0] = V1->Val[0]; 
      R->Val[4] = V1->Val[4]; 
      R->Val[8] = V1->Val[8];
      a1[0] = V1->Val[1];
      a1[1] = V1->Val[2];
      a1[2] = V1->Val[5];
      R->Val[1] = V1->Val[3]; 
      R->Val[2] = V1->Val[6];
      R->Val[5] = V1->Val[7];
      R->Val[3] = a1[0];
      R->Val[6] = a1[1];
      R->Val[7] = a1[2];
    }
    else{
      for(k=0 ; k<Current.NbrHar ; k+=1){
	R->Val[MAX_DIM*k+0] = V1->Val[MAX_DIM*k+0]; 
	R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4]; 
	R->Val[MAX_DIM*k+8] = V1->Val[MAX_DIM*k+8];
	a1[0] = V1->Val[MAX_DIM*k+1];
	a1[1] = V1->Val[MAX_DIM*k+2];
	a1[2] = V1->Val[MAX_DIM*k+5];
	R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+3]; 
	R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+6];
	R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+7];
	R->Val[MAX_DIM*k+3] = a1[0];
	R->Val[MAX_DIM*k+6] = a1[1];
	R->Val[MAX_DIM*k+7] = a1[2];
      }
    }
    break;    
    
  default:
    Msg(GERROR, "Wrong argument in 'Cal_TransposeValue'");
    break;
  }  
  
  GetDP_End ;
}

/* ------------------------------------------------------------------------ 
   R <- Trace(V1) 
   ------------------------------------------------------------------------ */

void Cal_TraceValue(struct Value *V1, struct Value *R){
  int     k;

  GetDP_Begin("Cal_TraceValue");

  switch(V1->Type){

  case TENSOR_DIAG :
    if (Current.NbrHar == 1) {
      R->Val[0] = V1->Val[0]+V1->Val[1]+V1->Val[2];
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k  ]+
	                    V1->Val[MAX_DIM*k+1]+
	                    V1->Val[MAX_DIM*k+2];
      }
    }
    R->Type = SCALAR ;
    break;

  case TENSOR_SYM :
    if (Current.NbrHar == 1) {
      R->Val[0] = V1->Val[0]+V1->Val[3]+V1->Val[5];
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k  ]+
	                    V1->Val[MAX_DIM*k+3]+
	                    V1->Val[MAX_DIM*k+5];
      }
    }
    R->Type = SCALAR ;
    break;

  case TENSOR :
    if (Current.NbrHar == 1) {
      R->Val[0] = V1->Val[0]+V1->Val[4]+V1->Val[8];
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k  ]+
	                    V1->Val[MAX_DIM*k+4]+
	                    V1->Val[MAX_DIM*k+8];
      }
    }
    R->Type = SCALAR ;
    break;
    
  default:
    Msg(GERROR, "Wrong argument type in 'Cal_TraceValue'");
    break;
  }  

  GetDP_End ;  
}

/* ------------------------------------------------------------------------ 
   R <- V1^T * V2 * V1  ,  V1 real 
   ------------------------------------------------------------------------ */

#define A0  V1->Val[0]
#define A1  V1->Val[1]
#define A2  V1->Val[2]
#define A3  V1->Val[3]
#define A4  V1->Val[4]
#define A5  V1->Val[5]
#define A6  V1->Val[6]
#define A7  V1->Val[7]
#define A8  V1->Val[8]

void Cal_RotateValue(struct Value *V1, struct Value *V2, struct Value *R){
  int           k;
  double        t0,t1,t2,t3,t4,t5,t6,t7,t8;
  struct Value  A;

  GetDP_Begin("Cal_RotateValue");

  switch(V2->Type){

  case VECTOR :
    if(Current.NbrHar == 1){
#define B0  V2->Val[0]
#define B1  V2->Val[1]
#define B2  V2->Val[2]
      A.Val[0]= A0*B0+A1*B1+A2*B2;
      A.Val[1]= A3*B0+A4*B1+A5*B2;
      A.Val[2]= A6*B0+A7*B1+A8*B2;
      A.Type = VECTOR;
      Cal_CopyValue(&A,R); 
#undef B0
#undef B1
#undef B2
    }
    else{ /* Attention: a modifier */
#define B0  V2->Val[0]
#define B1  V2->Val[1]
#define B2  V2->Val[2]
      A.Val[0]= A0*B0+A1*B1+A2*B2;
      A.Val[1]= A3*B0+A4*B1+A5*B2;
      A.Val[2]= A6*B0+A7*B1+A8*B2;
      A.Type = VECTOR;
      Cal_CopyValue(&A,R); 
#undef B0
#undef B1
#undef B2
    }
    break ;

  case TENSOR_DIAG :
    if(Current.NbrHar == 1){
#define B0  V2->Val[0]
#define B1  V2->Val[1]
#define B2  V2->Val[2]
      A.Val[0]= A0*A0*B0+A3*A3*B1+A6*A6*B2;
      A.Val[1]= A1*A0*B0+A3*A4*B1+A6*A7*B2;
      A.Val[2]= A2*A0*B0+A3*A5*B1+A6*A8*B2;
      A.Val[3]= A1*A0*B0+A3*A4*B1+A6*A7*B2;
      A.Val[4]= A1*A1*B0+A4*A4*B1+A7*A7*B2;
      A.Val[5]= A2*A1*B0+A4*A5*B1+A7*A8*B2;
      A.Val[6]= A2*A0*B0+A3*A5*B1+A6*A8*B2;
      A.Val[7]= A2*A1*B0+A4*A5*B1+A7*A8*B2;
      A.Val[8]= A2*A2*B0+A5*A5*B1+A8*A8*B2;
      A.Type = TENSOR;
      Cal_CopyValue(&A,R); 
#undef B0
#undef B1
#undef B2
    }
    else{
#define B0r  V2->Val[MAX_DIM* k   +0]
#define B1r  V2->Val[MAX_DIM* k   +1]
#define B2r  V2->Val[MAX_DIM* k   +2]
#define B0i  V2->Val[MAX_DIM*(k+1)+0]
#define B1i  V2->Val[MAX_DIM*(k+1)+1]
#define B2i  V2->Val[MAX_DIM*(k+1)+2]
#define AFFECT(i)				\
  A.Val[MAX_DIM* k   +i] = t0*B0r+t1*B1r+t2*B2r;	\
  A.Val[MAX_DIM*(k+1)+i] = t0*B0i+t1*B1i+t2*B2i
      for(k=0 ; k<Current.NbrHar ; k+=2){	
	t0=A0*A0; t1=A3*A3; t2=A6*A6; AFFECT(0);
	t0=A1*A0; t1=A3*A4; t2=A6*A7; AFFECT(1);
	t0=A2*A0; t1=A3*A5; t2=A6*A8; AFFECT(2);
	t0=A1*A0; t1=A3*A4; t2=A6*A7; AFFECT(3);
	t0=A1*A1; t1=A4*A4; t2=A7*A7; AFFECT(4);
	t0=A2*A1; t1=A4*A5; t2=A7*A8; AFFECT(5);
	t0=A2*A0; t1=A3*A5; t2=A6*A8; AFFECT(6);
	t0=A2*A1; t1=A4*A5; t2=A7*A8; AFFECT(7);
	t0=A2*A2; t1=A5*A5; t2=A8*A8; AFFECT(8);	
      }
      A.Type = TENSOR;
      Cal_CopyValue(&A,R); 
#undef B0r
#undef B1r
#undef B2r
#undef B0i
#undef B1i
#undef B2i
#undef AFFECT
    }    
    break;

#define COMPUTE_A										\
  A.Val[0]= A0*A0*B0+A0*A3*B3+A0*A6*B6+A3*A0*B1+A3*A3*B4+A3*A6*B7+A6*A0*B2+A6*A3*B5+A6*A6*B8;	\
  A.Val[1]= A1*A0*B0+A1*A3*B3+A1*A6*B6+A4*A0*B1+A4*A3*B4+A4*A6*B7+A7*A0*B2+A7*A3*B5+A7*A6*B8;	\
  A.Val[2]= A2*A0*B0+A2*A3*B3+A2*A6*B6+A5*A0*B1+A5*A3*B4+A5*A6*B7+A8*A0*B2+A8*A3*B5+A8*A6*B8;	\
  A.Val[3]= A1*A0*B0+A0*A4*B3+A0*A7*B6+A3*A1*B1+A4*A3*B4+A3*A7*B7+A6*A1*B2+A6*A4*B5+A7*A6*B8;	\
  A.Val[4]= A1*A1*B0+A1*A4*B3+A1*A7*B6+A4*A1*B1+A4*A4*B4+A4*A7*B7+A7*A1*B2+A7*A4*B5+A7*A7*B8;	\
  A.Val[5]= A2*A1*B0+A2*A4*B3+A2*A7*B6+A5*A1*B1+A5*A4*B4+A5*A7*B7+A8*A1*B2+A8*A4*B5+A8*A7*B8;	\
  A.Val[6]= A2*A0*B0+A0*A5*B3+A0*A8*B6+A3*A2*B1+A5*A3*B4+A3*A8*B7+A6*A2*B2+A6*A5*B5+A8*A6*B8;	\
  A.Val[7]= A2*A1*B0+A1*A5*B3+A1*A8*B6+A4*A2*B1+A5*A4*B4+A4*A8*B7+A7*A2*B2+A7*A5*B5+A8*A7*B8;	\
  A.Val[8]= A2*A2*B0+A2*A5*B3+A2*A8*B6+A5*A2*B1+A5*A5*B4+A5*A8*B7+A8*A2*B2+A8*A5*B5+A8*A8*B8

#define COMPLE