pos_quantity.c

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

    Cal_PostQuantity(PostQuantity_P, DefineQuantity_P0, QuantityStorage_P0,
		     Support_L, &Element, 0, 0, 0, &(*Values)[i]) ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  C o m b i n e _ P o s t Q u a n t i t y                                 */
/* ------------------------------------------------------------------------ */

void Combine_PostQuantity(int Type, int Order, 
			  struct Value *V1, struct Value *V2){
  
  GetDP_Begin("Combine_PostQuantity");
  
  switch(Type){
  case MULTIPLICATION : Cal_ProductValue(V1, V2, V1) ; break ;
  case ADDITION :       Cal_AddValue(V1, V2, V1) ; break ;
  case DIVISION :       Cal_DivideValue(Order?V1:V2, Order?V2:V1, V1) ; break;
  case SOUSTRACTION :   Cal_SubstractValue(Order?V1:V2, Order?V2:V1, V1) ; break;	
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  P o s _ L o c a l O r I n t e g r a l Q u a n t i t y                   */
/* ------------------------------------------------------------------------ */

static int Warning_NoJacobian = 0 ;

void Pos_LocalOrIntegralQuantity(struct PostQuantity    *PostQuantity_P, 
				 struct DefineQuantity  *DefineQuantity_P0,
				 struct QuantityStorage *QuantityStorage_P0,
				 struct PostQuantityTerm  *PostQuantityTerm_P,
				 struct Element         *Element, 
				 int    Type_Quantity,
				 double u, double v, double w, 
				 struct Value *Value) {

  struct FunctionSpace     *FunctionSpace_P ;
  struct DefineQuantity    *DefineQuantity_P ;
  struct QuantityStorage   *QuantityStorage_P ;
  struct Value              TermValue, tmpValue;
  struct JacobianCase      *JacobianCase_P0 ;
  struct IntegrationCase   *IntegrationCase_P ;
  struct Quadrature        *Quadrature_P ;

  List_T   *IntegrationCase_L, *JacobianCase_L ;
  double    ui, vi, wi, weight, Factor ;
  int       Index_DefineQuantity ;
  int       i, j, Type_Dimension ;
  int       CriterionIndex, Nbr_IntPoints, i_IntPoint ;

  double   (*Get_Jacobian) (struct Element * Element, MATRIX3x3 * Jac) = 0;
  void     (*Get_IntPoint) (int Nbr_Points, int Num,
			    double * u, double * v, double * w, double * wght) ;

  GetDP_Begin("Pos_LocalOrIntegralQuantity");

  /* Get the functionspaces
     Get the DoF for local quantities */

  for (i = 0 ; i < PostQuantityTerm_P->NbrQuantityIndex ; i++) {
    Index_DefineQuantity = PostQuantityTerm_P->QuantityIndexTable[i] ;
    DefineQuantity_P     = DefineQuantity_P0  + Index_DefineQuantity ;
    QuantityStorage_P    = QuantityStorage_P0 + Index_DefineQuantity ;
	
    if (QuantityStorage_P->NumLastElementForFunctionSpace != Element->Num) {

      QuantityStorage_P->NumLastElementForFunctionSpace = Element->Num ;

      if (Type_Quantity != INTEGRALQUANTITY){
	QuantityStorage_P->FunctionSpace = FunctionSpace_P = 
	  (struct FunctionSpace*)
	  List_Pointer(Problem_S.FunctionSpace,
		       DefineQuantity_P->FunctionSpaceIndex) ;
	if (DefineQuantity_P->Type == LOCALQUANTITY)
	  Get_DofOfElement(Element, FunctionSpace_P, QuantityStorage_P,
			   DefineQuantity_P->IndexInFunctionSpace) ;
      }
      else{ /* INTEGRALQUANTITY */
	if(DefineQuantity_P->IntegralQuantity.DefineQuantityIndexDof >= 0)
	  QuantityStorage_P->FunctionSpace = (struct FunctionSpace*)
	    List_Pointer(Problem_S.FunctionSpace,
			 DefineQuantity_P->FunctionSpaceIndex) ;
	/* Get the function space for the associated local quantities */
	for (j = 0 ; j < DefineQuantity_P->IntegralQuantity.NbrQuantityIndex ; j++) {
	  Index_DefineQuantity = DefineQuantity_P->IntegralQuantity.QuantityIndexTable[j];
	  DefineQuantity_P     = DefineQuantity_P0  + Index_DefineQuantity ;
	  QuantityStorage_P    = QuantityStorage_P0 + Index_DefineQuantity ;
	  QuantityStorage_P->FunctionSpace = (struct FunctionSpace*)
	    List_Pointer(Problem_S.FunctionSpace,
			 DefineQuantity_P->FunctionSpaceIndex) ;	    
	}
      }

    }
  }

  /* get the jacobian */

  if (Element->Num != NO_ELEMENT) {
    if (PostQuantityTerm_P->JacobianMethodIndex >= 0) {	  
      JacobianCase_L = 
	((struct JacobianMethod *)
	 List_Pointer(Problem_S.JacobianMethod, 
		      PostQuantityTerm_P->JacobianMethodIndex))
	->JacobianCase ;
      JacobianCase_P0 = (struct JacobianCase*)List_Pointer(JacobianCase_L, 0);

      i = 0 ;
      while ((i < List_Nbr(JacobianCase_L)) &&
	     ((j = (JacobianCase_P0 + i)->RegionIndex) >= 0) &&
	     !List_Search
	     (((struct Group *)List_Pointer(Problem_S.Group, j))
	      ->InitialList, &Element->Region, fcmp_int) )  i++ ;

      if (i == List_Nbr(JacobianCase_L))
	Msg(GERROR, "Undefined Jacobian in Region %d", Element->Region) ;

      Element->JacobianCase = JacobianCase_P0 + i ;
      Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
	Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
			     Element->Type, &Type_Dimension) ;
    }
    else {
      if(!Warning_NoJacobian){
	Msg(WARNING, "No Jacobian method specification in PostProcessing quantity");
	Msg(WARNING, "Using default Jacobian (Vol)");
	Warning_NoJacobian = 1 ;
      }
      Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
	Get_JacobianFunction(JACOBIAN_VOL, Element->Type, &Type_Dimension) ;
    }

    Get_NodesCoordinatesOfElement(Element) ;
  }

  /* local evaluation at one point */

  if(PostQuantityTerm_P->EvaluationType == LOCAL){

    if (Element->Num != NO_ELEMENT) {
      Get_BFGeoElement(Element, u, v, w) ;
      Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
      if (Element->DetJac != 0.)
	Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
			  &Element->Jac, &Element->InvJac) ;	  
    }

    Cal_WholeQuantity
      (Current.Element = Element,
       QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
       Current.u = u, Current.v = v, Current.w = w, -1, -1, &TermValue) ;

  }

  /* integral evaluation over the element */

  else if(PostQuantityTerm_P->EvaluationType == INTEGRAL){

    if(Element->Num == NO_ELEMENT)
      Msg(GERROR, "No element in which to integrate");

    if(PostQuantityTerm_P->IntegrationMethodIndex < 0)
      Msg(GERROR, "Missing Integration method in PostProcesssing Quantity '%s'", 
	  PostQuantity_P->Name);
    
    IntegrationCase_L = 
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod, 
		    PostQuantityTerm_P->IntegrationMethodIndex))->IntegrationCase ;
    
    CriterionIndex = 
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod, 
		    PostQuantityTerm_P->IntegrationMethodIndex))
      ->CriterionIndex ;
    
    IntegrationCase_P = Get_IntegrationCase(Element, 
					    IntegrationCase_L,
					    CriterionIndex) ;
    
    if(IntegrationCase_P->Type != GAUSS)
      Msg(GERROR, "Only numerical integration is available "
	  "in Integral PostQuantities");
    
    Quadrature_P = (struct Quadrature*)
      List_PQuery(IntegrationCase_P->Case, &Element->Type, fcmp_int);
    
    if(!Quadrature_P)
      Msg(GERROR, "Unknown type of Element (%s) for Integration method (%s) "
	  " in PostProcessing Quantity (%s)", 
	  Get_StringForDefine(Element_Type, Element->Type),
	  ((struct IntegrationMethod *)
	   List_Pointer(Problem_S.IntegrationMethod,
			PostQuantityTerm_P->IntegrationMethodIndex))->Name,
	  PostQuantity_P->Name);
    
    Cal_ZeroValue(&TermValue);
    
    Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
    Get_IntPoint = (void (*) (int,int,double*,double*,double*,double*))
      Quadrature_P->Function ;
    
    for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
      
      Get_IntPoint(Nbr_IntPoints, i_IntPoint, &ui, &vi, &wi, &weight) ;
      
      Get_BFGeoElement (Element, ui, vi, wi) ;
      Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
      if (Element->DetJac != 0.){
	Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
			  &Element->Jac, &Element->InvJac) ;	  
      }
      else{
	Msg(WARNING, "Zero determinant in 'Cal_PostQuantity'");
      }
      Current.x = Current.y = Current.z = 0. ;
      if (Type_Quantity == INTEGRALQUANTITY){
	for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
	  Current.x += Element->x[i] * Element->n[i] ;
	  Current.y += Element->y[i] * Element->n[i] ;
	  Current.z += Element->z[i] * Element->n[i] ;
	}
      }
      
      Cal_WholeQuantity
	(Current.Element = Element,
	 QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
	 Current.u = ui, Current.v = vi, Current.w = wi, -1, -1, &tmpValue) ;
      
      Factor = weight * fabs(Element->DetJac) ;
      
      TermValue.Type = tmpValue.Type ;
      Cal_AddMultValue(&TermValue,&tmpValue,Factor,&TermValue);
    }
  }

  Value->Type = TermValue.Type;
  Cal_AddValue(Value,&TermValue,Value);

  GetDP_End ;
}