cal_quantity.c

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

	xyz2uvwInAnElement(Element->ElementTrace, Current.x, Current.y, Current.z, 
			   &Current.ut, &Current.vt, &Current.wt) ;	
	Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0,
			  WholeQuantity_P->Case.Trace.WholeQuantity,
			  Current.ut, Current.vt, Current.wt, 
			  -1, 0, &Stack[0][Index]) ;
      }

      Current.Region = Save_Region ;
      Multi[Index] = 0 ;
      Index++ ;
      break ;

    case WQ_SOLIDANGLE : /* SolidAngle[{qty}] */
      Cal_SolidAngle(0, Element, QuantityStorage_P0 +
		     WholeQuantity_P->Case.OperatorAndQuantity.Index, 
		     Nbr_Dof, Index, (struct Value **)Stack);
      Multi[Index] = 1 ;
      Index++ ;
      break ;

    case WQ_BINARYOPERATOR : /* + - * x / % ^ < > <= >= == != && || */
      if (Index-2 != DofIndex && Index-1 != DofIndex){
	if(!Multi[Index-1] && !Multi[Index-2])
	  ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	    (&Stack[0][Index-2], &Stack[0][Index-1], &Stack[0][Index-2]) ;
	else if(Multi[Index-1] && Multi[Index-2])
	  for(j=0 ; j<Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&Stack[j][Index-2], &Stack[j][Index-1], &Stack[j][Index-2]) ;
	else if(Multi[Index-2])
	  for(j=0 ; j<Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&Stack[j][Index-2], &Stack[0][Index-1], &Stack[j][Index-2]) ;
	else {
	  for(j=0 ; j<Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&Stack[0][Index-2], &Stack[j][Index-1], &Stack[j][Index-2]) ;
	  Multi[Index-2] = 1 ;
	}
      }
      else if (Index-1 == DofIndex) {
	if(Multi[Index-2])
	  for (j = 0 ; j < Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&Stack[j][Index-2], &DofValue[j], &DofValue[j]) ;
	else
	  for (j = 0 ; j < Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&Stack[0][Index-2], &DofValue[j], &DofValue[j]) ;
	DofIndex-- ;
      }
      else { /* Index-2 == DofIndex */
	if(Multi[Index-1])
	  for (j = 0 ; j < Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&DofValue[j], &Stack[j][Index-1], &DofValue[j]) ;
	else
	  for (j = 0 ; j < Nbr_Dof ; j++)
	    ((CAST3V)WholeQuantity_P->Case.Operator.Function)
	      (&DofValue[j], &Stack[0][Index-1], &DofValue[j]) ;
      }
      Index-- ;  
      break ;

    case WQ_UNARYOPERATOR : /* + - ! */
      if (Index-1 == DofIndex)
	for (j = 0 ; j < Nbr_Dof ; j++)
	  ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&DofValue[j]) ;
      else if(Multi[Index-1])
	for(j=0 ; j<Nbr_Dof ; j++)
	  ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[j][Index-1]) ;	
      else
	((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[0][Index-1]) ;
      break ;

      /* WARNING: all the rest assumes 0 multi status */

    case WQ_TEST :
      Flag_True = (Stack[0][Index-1].Val[0] != 0.) ;
      Cal_WholeQuantity(Element, QuantityStorage_P0,
			(Flag_True) ? WholeQuantity_P->Case.Test.WholeQuantity_True :
		 	              WholeQuantity_P->Case.Test.WholeQuantity_False,
			u, v, w, -1, 0, &Stack[0][Index-1]) ;
      break ;

    case WQ_EXPRESSION :
      Index -= WholeQuantity_P->Case.Expression.NbrArguments ;
      Get_ValueOfExpression
	(Problem_Expression0 + WholeQuantity_P->Case.Expression.Index, 
	 QuantityStorage_P0, u, v, w, &Stack[0][Index]) ;
      Multi[Index] = 0 ;
      Index++ ;
      break ;

    case WQ_BUILTINFUNCTION :
      Index -= WholeQuantity_P->Case.Function.NbrArguments ;
      ((CASTF2V)WholeQuantity_P->Case.Function.Fct)
	(&WholeQuantity_P->Case.Function, &Stack[0][Index], &Stack[0][Index]) ;	
      Multi[Index] = 0 ;
      Index++ ;
      break ;

    case WQ_EXTERNBUILTINFUNCTION :
      ((CASTF2V)WholeQuantity_P->Case.Function.Fct) 
	(&WholeQuantity_P->Case.Function, DofValue, &Stack[0][Index]) ;
      Multi[Index] = 0 ;
      Index++ ;
      break ;

    case WQ_CONSTANT :
      if (Current.NbrHar == 1) {
	Stack[0][Index].Val[0] = WholeQuantity_P->Case.Constant ;
      }
      else {
	for (k = 0 ; k < Current.NbrHar ; k += 2) {
	  Stack[0][Index].Val[MAX_DIM* k   ] = WholeQuantity_P->Case.Constant ;
	  Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ;
	}
      }
      Stack[0][Index].Type = SCALAR ;
      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_CURRENTVALUE :
      if (Current.NbrHar == 1) {
	Stack[0][Index].Val[0] = *(WholeQuantity_P->Case.CurrentValue.Value) ;
      }
      else {
	for (k = 0 ; k < Current.NbrHar ; k += 2) {
	  Stack[0][Index].Val[MAX_DIM* k   ] = *(WholeQuantity_P->Case.CurrentValue.Value) ;
	  Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ;
	}
      }
      Stack[0][Index].Type = SCALAR ;
      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_ARGUMENT :
      Cal_CopyValue(DofValue + WholeQuantity_P->Case.Argument.Index - 1,
		    &Stack[0][Index]) ;
      Multi[Index] = 0 ;
      Index++ ;
      break ;

    case WQ_TIMEDERIVATIVE :
      if (Current.NbrHar == 1) {
	Cal_WholeQuantity(Element, QuantityStorage_P0,
			  WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
			  u, v, w, -1, 0, &Stack[0][Index]) ;
	
	for (k = 0 ; k < Current.NbrSystem ; k++){
	  if(List_Nbr((Current.DofData_P0+k)->Solutions) > 1){
	    Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)->Solutions, 0);
	    if ((Current.DofData_P0+k)->CurrentSolution != Solution_P0)
	      ((Current.DofData_P0+k)->CurrentSolution) -- ;
	  }
	  else {
	    if (!Flag_WarningMissSolForDt) {
	      Msg(WARNING,
		  "Missing solution for time derivative computation (Sys#%d/%d)",
		  k, Current.NbrSystem);
	      Flag_WarningMissSolForDt = 1 ;
	    }
	  }
	}

	Save_Time = Current.Time ;
	Current.Time = Current.DofData->CurrentSolution->Time ;

	Cal_WholeQuantity(Element, QuantityStorage_P0,
			  WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
			  u, v, w, -1, 0, &Stack[0][Index+1]) ;
	Cal_SubstractValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) ;
	Stack[0][Index+1].Val[0] = Save_Time - Current.Time ;
	Stack[0][Index+1].Type = SCALAR ;
	if(Stack[0][Index+1].Val[0])
	  Cal_DivideValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) ;
	else
	  Cal_ZeroValue(&Stack[0][Index]);

	for (k = 0 ; k < Current.NbrSystem ; k++)
	  if (List_Nbr((Current.DofData_P0+k)->Solutions) > 1)
	    ((Current.DofData_P0+k)->CurrentSolution) ++ ;

	Current.Time = Save_Time ;
      }
      else {
	Cal_WholeQuantity(Element, QuantityStorage_P0,
			  WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
			  u, v, w, -1, 0, &Stack[0][Index]) ;
	for (k = 0 ; k < Current.NbrHar ; k += 2) {
	  Stack[0][Index+1].Val[MAX_DIM* k   ] = 0. ;
	  Stack[0][Index+1].Val[MAX_DIM*(k+1)] = Current.DofData->Val_Pulsation[k/2] ;
	}
	Stack[0][Index+1].Type = SCALAR ;
	Cal_ProductValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) ;
      }
      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_ATANTERIORTIMESTEP :
      ntime = WholeQuantity_P->Case.AtAnteriorTimeStep.TimeStep ;
      
      for (k = 0 ; k < Current.NbrSystem ; k++){
	Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)->Solutions, 0);
	if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) >= ntime){ 
	  ((Current.DofData_P0+k)->CurrentSolution) -= ntime ;
	  if (Flag_InfoForTime_ntime != List_Nbr((Current.DofData_P0+k)->Solutions)) {
	    Msg(INFO, "Accessing solution from %d time steps ago", ntime);
	    Msg(INFO, "  -> System %d/%d: TimeStep = %d, Time = %g + i * %g",
		k+1, Current.NbrSystem, 
		(Current.DofData_P0+k)->CurrentSolution->TimeStep,
		(Current.DofData_P0+k)->CurrentSolution->Time,
		(Current.DofData_P0+k)->CurrentSolution->TimeImag);
	    Flag_InfoForTime_ntime = List_Nbr((Current.DofData_P0+k)->Solutions);
	  }
	}
	else {
	  if (!Flag_WarningMissSolForTime_ntime) {
	    Msg(WARNING,
		"Missing solution for time step -%d computation (System #%d/%d)",
		ntime, k+1, Current.NbrSystem);
	    Flag_WarningMissSolForTime_ntime = 1 ;
	  }
	}
      }
      
      Save_TimeStep = Current.TimeStep ;
      Save_Time = Current.Time ;
      Save_TimeImag = Current.TimeImag ;
      Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ;
      Current.Time = Current.DofData->CurrentSolution->Time ;
      Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ;
      
      Cal_WholeQuantity(Element, QuantityStorage_P0,
			WholeQuantity_P->Case.AtAnteriorTimeStep.WholeQuantity,
			u, v, w, -1, 0, &Stack[0][Index]) ;
      
      Current.TimeStep = Save_TimeStep ;
      Current.Time = Save_Time ;
      Current.TimeImag = Save_TimeImag ;

      for (k = 0 ; k < Current.NbrSystem ; k++){
	Solution_PN = (struct Solution*)
	  List_Pointer((Current.DofData_P0+k)->Solutions, 
		       List_Nbr((Current.DofData_P0+k)->Solutions)-1);
	if((Solution_PN - (Current.DofData_P0+k)->CurrentSolution) >= ntime)
	  ((Current.DofData_P0+k)->CurrentSolution) += ntime ;
      }

      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_MHTRANSFORM :
      if(Current.NbrHar == 1)
	Msg(GERROR, 
	    "MHTransform can only be used in complex (multi-harmonic) calculations") ;

      Cal_WholeQuantity(Element, QuantityStorage_P0, 
			WholeQuantity_P->Case.MHTransform.WholeQuantity, 
			u, v, w, -1, 0, &Stack[0][Index]) ; 

      MHTransform(Element, QuantityStorage_P0, u, v, w, &Stack[0][Index],
		  Problem_Expression0 + WholeQuantity_P->Case.MHTransform.Index,
		  WholeQuantity_P->Case.MHTransform.NbrPoints) ;
      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_CAST : 
      /* This should be changed... */
      Save_NbrHar = Current.NbrHar ;
      Save_DofData = Current.DofData ;
         
      if (!WholeQuantity_P->Case.Cast.NbrHar){
	Current.DofData = 	
	  ((struct FunctionSpace *)
	   List_Pointer(Problem_S.FunctionSpace,
			WholeQuantity_P->Case.Cast.FunctionSpaceIndexForType))
	  ->DofData ;
	
	Current.NbrHar = Current.DofData->NbrHar ;
      }
      else{
	Current.NbrHar = WholeQuantity_P->Case.Cast.NbrHar ;
      }

      Cal_WholeQuantity(Element, QuantityStorage_P0,
			WholeQuantity_P->Case.Cast.WholeQuantity,
			u, v, w, -1, 0, &Stack[0][Index]) ;
      
      if (Current.NbrHar < Save_NbrHar)  /* ne plus a completer ...?? */
      	Cal_SetZeroHarmonicValue(&Stack[0][Index], Save_NbrHar) ;
      
      Current.NbrHar = Save_NbrHar ;
      Current.DofData = Save_DofData ;
      Multi[Index] = 0 ;
      Index++ ;
      break ;

    case WQ_SAVEVALUE :
      if(WholeQuantity_P->Case.SaveValue.Index > MAX_REGISTER_SIZE-1)
	Msg(GERROR, "Register Size Exceeded (%d)", MAX_REGISTER_SIZE);
      /* if (WholeQuantity_P->Case.SaveValue.Index >= 0) */
      Cal_CopyValue(&Stack[0][Index-1], 
		    ValueSaved + WholeQuantity_P->Case.SaveValue.Index) ;
      break ;

    case WQ_VALUESAVED :
      if(WholeQuantity_P->Case.ValueSaved.Index > MAX_REGISTER_SIZE-1)
	Msg(GERROR, "Register size exceeded (%d)", MAX_REGISTER_SIZE);
      Cal_CopyValue(ValueSaved + WholeQuantity_P->Case.ValueSaved.Index, 
		    &Stack[0][Index]) ;
      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_SHOWVALUE :
      if (Index-1 == DofIndex) {
	for(j=0 ; j<Nbr_Dof ; j++){
	  fprintf(stderr, "##%d Dof %d ", WholeQuantity_P->Case.ShowValue.Index, j+1);
	  Show_Value(&DofValue[j]);
	}
      } 
      else {
	fprintf(stderr, "##%d ", WholeQuantity_P->Case.ShowValue.Index);
	Show_Value(&Stack[0][Index-1]);
      }
      break ;
      
    default :
      Msg(GERROR, "Unknown type of WholeQuantity (%d)", WholeQuantity_P->Type);
      break;
    }
  }

  if (DofIndexInWholeQuantity < 0) Cal_CopyValue(&Stack[0][0], &DofValue[0]) ;

#if defined(USE_GLOBAL_STACK)
  RecursionIndex--;
#endif

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  P u r i f y _ W h o l e Q u a n t i t y                                 */
/* ------------------------------------------------------------------------ */

List_T * Purify_WholeQuantity(List_T * WQ_L) {
  
  GetDP_Begin("Purify_WholeQuantity");

  /* It would be nice to pre-compute all trivial sequences in a list
     of WholeQuantities. For example, when all the quantties are
     constants, it is pretty stupid to recompute everything
     everytime using the stack... */

  GetDP_Return(NULL) ;
}

/* ------------------------------------------------------------------------ */
/*  C a l _ S t o r e I n R e g i s t e r                                   */
/* ------------------------------------------------------------------------ */

void Cal_StoreInRegister(struct  Value  *Value, int RegisterIndex ) {
  Cal_CopyValue(Value, ValueSaved + RegisterIndex) ;
}