cal_value.c

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

	R[i].Type = V1[i].Type;
	for (k = 0 ; k < Current.NbrHar ; k++) {
	  CADD(i,0); CADD(i,1); CADD(i,2);
	}
      }
  }
  else if (V1[0].Type == TENSOR_SYM && V2[0].Type == TENSOR_SYM){
    if (Current.NbrHar == 1)
      for(i = 0 ; i < Nbr_Values; i++){ 
	R[i].Type = TENSOR_SYM;
	ADD(i,0); ADD(i,1); ADD(i,2); ADD(i,3); ADD(i,4); ADD(i,5);
      }
    else
      for(i = 0 ; i < Nbr_Values; i++){ 
	R[i].Type = TENSOR_SYM;
	for (k = 0 ; k < Current.NbrHar ; k++) {
	  CADD(i,0); CADD(i,1); CADD(i,2); CADD(i,3); CADD(i,4); CADD(i,5);
	}
      }
  }
  else if (V1[0].Type == TENSOR && V2[0].Type == TENSOR){
    if (Current.NbrHar == 1)
      for(i = 0 ; i < Nbr_Values; i++){
	R[i].Type = TENSOR;
	ADD(i,0); ADD(i,1); ADD(i,2); ADD(i,3); ADD(i,4); ADD(i,5); ADD(i,6); ADD(i,7); ADD(i,8);
      }
    else
      for(i = 0 ; i < Nbr_Values; i++){
	R[i].Type = TENSOR;
	for (k = 0 ; k < Current.NbrHar ; k++) {
	  CADD(i,0); CADD(i,1); CADD(i,2); CADD(i,3); CADD(i,4); CADD(i,5); CADD(i,6); CADD(i,7); CADD(i,8);
	}
      }
  }
  else if ((V1[0].Type == TENSOR     && V2[0].Type == TENSOR_SYM) ||
	   (V1[0].Type == TENSOR     && V2[0].Type == TENSOR_DIAG)||
	   (V1[0].Type == TENSOR_SYM && V2[0].Type == TENSOR_DIAG)){
    A.Type = V1[0].Type;
    for(i = 0 ; i < Nbr_Values; i++){
      for (k = 0 ; k < Current.NbrHar ; k++) {
	for(ii=0 ; ii<9 ; ii++){
	  i1 = (V1[0].Type==TENSOR)?ii:TENSOR_SYM_MAP[ii];
	  i2 = (V2[0].Type==TENSOR_SYM)?TENSOR_SYM_MAP[ii]:TENSOR_DIAG_MAP[ii];
	  A.Val[MAX_DIM*k+i1] = 
	    V1[i].Val[MAX_DIM*k+i1] + ((i2<0)?0.0:V2[i].Val[MAX_DIM*k+i2]);
	}
      }
      Cal_CopyValue(&A,&R[i]);
    }
  }
  else if ((V1[0].Type == TENSOR_SYM  && V2[0].Type == TENSOR)      ||
	   (V1[0].Type == TENSOR_DIAG && V2[0].Type == TENSOR)      ||
	   (V1[0].Type == TENSOR_DIAG && V2[0].Type == TENSOR_SYM)){
    A.Type = V2[0].Type;
    for(i = 0 ; i < Nbr_Values; i++){
      for (k = 0 ; k < Current.NbrHar ; k++) {
	for(ii=0 ; ii<9 ; ii++){
	  i1 = (V1[i].Type==TENSOR_SYM)?TENSOR_SYM_MAP[ii]:TENSOR_DIAG_MAP[ii];
	  i2 = (V2[i].Type==TENSOR)?ii:TENSOR_SYM_MAP[ii];
	  A.Val[MAX_DIM*k+i2] = 
	    ((i1<0)?0.0:V1[i].Val[MAX_DIM*k+i1]) + V2[i].Val[MAX_DIM*k+i2];
	}
      }
      Cal_CopyValue(&A,&R[i]);
    }
  }
  else
    Msg(GERROR, "Addition of different quantities: %s + %s",
	Get_StringForDefine(Field_Type, V1->Type),
	Get_StringForDefine(Field_Type, V2->Type));
  
  GetDP_End ;  
}

#undef ADD
#undef CADD

/* ------------------------------------------------------------------------ 
   R <- V1 * d ,   where d is a double
   ------------------------------------------------------------------------ */

void  Cal_MultValue (struct Value * V1, double d, struct Value * R) {
  int k;

  GetDP_Begin("Cal_MultValue");

  R->Type = V1->Type ;

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



/* ------------------------------------------------------------------------ 
   R <- V1 + V2 * d ,   where d is a double
   ------------------------------------------------------------------------ */

void  Cal_AddMultValue (struct Value * V1, struct Value * V2, double d, struct Value * R) {
  int k;
  struct Value A ;

  GetDP_Begin("Cal_AddMultValue");

  A.Type = V2->Type ;

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

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- V1 + V2 * d ,   where d is a double
   ------------------------------------------------------------------------ */

void  Cal_AddMultValueArray (struct Value * V1, struct Value * V2, double d, struct Value * R,int Nbr_Values) {
  int k,i;

  GetDP_Begin("Cal_AddMultValueArray");

  switch(V2[0].Type){
  case SCALAR :
    if (Current.NbrHar == 1) {
      for(i = 0 ; i < Nbr_Values ; i++)
	V2[i].Val[0] = V2[i].Val[0] * d;      
    }
    else{
      for (k = 0 ; k < Current.NbrHar ; k++) {
	for(i = 0 ; i < Nbr_Values ; i++)
	  V2[i].Val[MAX_DIM*k] = V2[i].Val[MAX_DIM*k] * d;
      }
    }
    break;
  case VECTOR :
  case TENSOR_DIAG :
    if (Current.NbrHar == 1) {
      for(i = 0 ; i < Nbr_Values ; i++){
	V2[i].Val[0] = V2[i].Val[0] * d;
	V2[i].Val[1] = V2[i].Val[1] * d;
	V2[i].Val[2] = V2[i].Val[2] * d;
      }
    }
    else{
      for (k = 0 ; k < Current.NbrHar ; k++) {
	for(i = 0 ; i < Nbr_Values ; i++){
	  V2[i].Val[MAX_DIM*k  ] = V2[i].Val[MAX_DIM*k  ] * d;
	  V2[i].Val[MAX_DIM*k+1] = V2[i].Val[MAX_DIM*k+1] * d;
	  V2[i].Val[MAX_DIM*k+2] = V2[i].Val[MAX_DIM*k+2] * d;
	}
      }
    }
    break;
  case TENSOR_SYM :
    if (Current.NbrHar == 1) {
      for(i = 0 ; i < Nbr_Values ; i++){
	V2[i].Val[0] = V2[i].Val[0] * d;
	V2[i].Val[1] = V2[i].Val[1] * d;
	V2[i].Val[2] = V2[i].Val[2] * d;
	V2[i].Val[3] = V2[i].Val[3] * d;
	V2[i].Val[4] = V2[i].Val[4] * d;
	V2[i].Val[5] = V2[i].Val[5] * d;
      }
    }
    else{
      for (k = 0 ; k < Current.NbrHar ; k++) {
	for(i = 0 ; i < Nbr_Values ; i++){
	  V2[i].Val[MAX_DIM*k  ] = V2[i].Val[MAX_DIM*k  ] * d;
	  V2[i].Val[MAX_DIM*k+1] = V2[i].Val[MAX_DIM*k+1] * d;
	  V2[i].Val[MAX_DIM*k+2] = V2[i].Val[MAX_DIM*k+2] * d;
	  V2[i].Val[MAX_DIM*k+3] = V2[i].Val[MAX_DIM*k+3] * d;
	  V2[i].Val[MAX_DIM*k+4] = V2[i].Val[MAX_DIM*k+4] * d;
	  V2[i].Val[MAX_DIM*k+5] = V2[i].Val[MAX_DIM*k+5] * d;
	}
      }
    }
    break;
  case TENSOR :
    if (Current.NbrHar == 1) {
      for(i = 0 ; i < Nbr_Values ; i++){
	V2[i].Val[0] = V2[i].Val[0] * d;
	V2[i].Val[1] = V2[i].Val[1] * d;
	V2[i].Val[2] = V2[i].Val[2] * d;
	V2[i].Val[3] = V2[i].Val[3] * d;
	V2[i].Val[4] = V2[i].Val[4] * d;
	V2[i].Val[5] = V2[i].Val[5] * d;
	V2[i].Val[6] = V2[i].Val[6] * d;
	V2[i].Val[7] = V2[i].Val[7] * d;
	V2[i].Val[8] = V2[i].Val[8] * d;
      }
    }
    else{
      for (k = 0 ; k < Current.NbrHar ; k++) {
	for(i = 0 ; i < Nbr_Values ; i++){
	  V2[i].Val[MAX_DIM*k  ] = V2[i].Val[MAX_DIM*k  ] * d;
	  V2[i].Val[MAX_DIM*k+1] = V2[i].Val[MAX_DIM*k+1] * d;
	  V2[i].Val[MAX_DIM*k+2] = V2[i].Val[MAX_DIM*k+2] * d;
	  V2[i].Val[MAX_DIM*k+3] = V2[i].Val[MAX_DIM*k+3] * d;
	  V2[i].Val[MAX_DIM*k+4] = V2[i].Val[MAX_DIM*k+4] * d;
	  V2[i].Val[MAX_DIM*k+5] = V2[i].Val[MAX_DIM*k+5] * d;
	  V2[i].Val[MAX_DIM*k+6] = V2[i].Val[MAX_DIM*k+6] * d;
	  V2[i].Val[MAX_DIM*k+7] = V2[i].Val[MAX_DIM*k+7] * d;
	  V2[i].Val[MAX_DIM*k+8] = V2[i].Val[MAX_DIM*k+8] * d;
	}
      }
    }
    break;        
  default :
    Msg(GERROR, "Wrong argument type for 'Cal_AddMultValueArray'");
    break;
  }
  Cal_AddValueArray(V1, V2, R,Nbr_Values);
  
  GetDP_End ;
}






/* ------------------------------------------------------------------------ 
   V1 <- V1 * d1 + V2 * d2 ,   where d1, d2 are doubles
   ------------------------------------------------------------------------ */

void  Cal_AddMultValue2 (struct Value * V1, double d1,
			 struct Value * V2, double d2) {
  int k;

  GetDP_Begin("Cal_AddMultValue");

  switch(V1->Type){
  case SCALAR :
    if (Current.NbrHar == 1) {
      V1->Val[0] = V1->Val[0] * d1 + V2->Val[0] * d2;      
    }
    else{
      for (k = 0 ; k < Current.NbrHar ; k++) {
	V1->Val[MAX_DIM*k] = V1->Val[MAX_DIM*k] * d1 + V2->Val[MAX_DIM*k] * d2;
      }
    }
    break;
  case VECTOR :
  case TENSOR_DIAG :
    if (Current.NbrHar == 1) {
      V1->Val[0] = V1->Val[0] * d1 + V2->Val[0] * d2;
      V1->Val[1] = V1->Val[1] * d1 + V2->Val[1] * d2;
      V1->Val[2] = V1->Val[2] * d1 + V2->Val[2] * d2;
    }
    else{
      for (k = 0 ; k < Current.NbrHar ; k++) {
	V1->Val[MAX_DIM*k  ] = V1->Val[MAX_DIM*k  ] * d1 + V2->Val[MAX_DIM*k  ] * d2;
	V1->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] * d1 + V2->Val[MAX_DIM*k+1] * d2;
	V1->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] * d1 + V2->Val[MAX_DIM*k+2] * d2;
      }
    }
    break;
  default :
    Msg(GERROR, "Wrong argument type for 'Cal_AddMultValue2'");
    break;
  }

  GetDP_End ;
}

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

#define SUB(i)   R->Val[i] = V1->Val[i] - V2->Val[i]
#define CSUB(i)  R->Val[MAX_DIM*k+i] = V1->Val[MAX_DIM*k+i] - V2->Val[MAX_DIM*k+i]

void  Cal_SubstractValue (struct Value * V1, struct Value * V2, struct Value * R) {
  int           i, k;
  int           i1, i2;
  struct Value  A;

  GetDP_Begin("Cal_SubstractValue");

  if (V1->Type == SCALAR && V2->Type == SCALAR) {