pos_feminterpolation.c

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

	Cal_InitIntegralQuantity (TheElement_P, &IQA, QuantityStorage_P);
	
      }
      else
	break ;
      
    }
    
    /* -----
       Local 
       ----- */

    if (Type_DefineQuantity == LOCALQUANTITY) {
      
      if (TheElement_P->Num != NO_ELEMENT) {
	for (j = 0 ; j < Nbr_Dof ; j++) {
	  xFunctionBF[j]
	    (TheElement_P, QuantityStorage_P->BasisFunction[j].NumEntityInElement+1,
	     u, v, w, vBFu[j]) ;
	  ((void (*)(struct Element*, double*, double*))
	   xChangeOfCoordinates) (TheElement_P, vBFu[j], vBFxDof[j].Val) ;
	  /*
   	  printf("j %d , Num %d, Type_Form %d  change %d "
	         "vBFu[j] %e %e %e  vBFx[j]  %e %e %e\n", 
		 j,  QuantityStorage_P->BasisFunction[j].NumEntityInElement+1,
		 Type_Form,(Flag_ChangeOfCoordinates && TheElement_P->Num != NO_ELEMENT),
		 vBFu[j][0],vBFu[j][1],vBFu[j][2],
		 vBFxDof[j].Val[0],vBFxDof[j].Val[1],vBFxDof[j].Val[2]);
	  */
	}
      }
      
      /* interpolate (vBFxDof is real-valued) */
      
      switch (Type_Form) {
	
      case FORM0 :  case FORM3 :  case FORM3P :  case SCALAR :
	if(Current.NbrHar==1){
	  for (j = 0 ; j < Nbr_Dof ; j++){
	    if (Type_Quantity == QUANTITY_BF)
	      Val_Dof = (QuantityStorage_P->BasisFunction[j].Dof->Entity ==
			 Current.SubRegion)? 1. : 0. ;
	    else
	      Dof_GetRealDofValue
		(QuantityStorage_P->FunctionSpace->DofData,
		 QuantityStorage_P->BasisFunction[j].Dof,
		 &Val_Dof) ;
	    Val[0] += vBFxDof[j].Val[0] * Val_Dof ;
	  }
	}
	else{	
	  for (j = 0 ; j < Nbr_Dof ; j++){
	    for (k = 0 ; k < Current.NbrHar ; k+=2) {
	      if (Type_Quantity == QUANTITY_BF) {
		Val_Dof_r = (QuantityStorage_P->BasisFunction[j].Dof->Entity ==
			     Current.SubRegion)? 1. : 0. ;
		Val_Dof_i = 0. ;
	      }
	      else {
		Dof_GetComplexDofValue
		  (QuantityStorage_P->FunctionSpace->DofData,
		   QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT,
		   &Val_Dof_r, &Val_Dof_i) ;
	      }
	      Val[MAX_DIM*k]     += vBFxDof[j].Val[0] * Val_Dof_r ;
	      Val[MAX_DIM*(k+1)] += vBFxDof[j].Val[0] * Val_Dof_i ;
	    }
	  }
	}
	*Type_Value = SCALAR ;
	break ;
	
      case FORM1  :  case FORM1P :  case FORM2  :  case FORM2P :
      case FORM1S :  case FORM2S :
      case VECTOR :  case VECTORP :
	if(Current.NbrHar==1){
	  for (j = 0 ; j < Nbr_Dof ; j++){
	    if (Type_Quantity == QUANTITY_BF)
	      Val_Dof = (QuantityStorage_P->BasisFunction[j].Dof->Entity ==
			 Current.SubRegion)? 1. : 0. ;
	    else
	      Dof_GetRealDofValue
		(QuantityStorage_P->FunctionSpace->DofData,
		 QuantityStorage_P->BasisFunction[j].Dof,
		 &Val_Dof) ;
	    Val[0] += vBFxDof[j].Val[0] * Val_Dof ;
	    Val[1] += vBFxDof[j].Val[1] * Val_Dof ;
	    Val[2] += vBFxDof[j].Val[2] * Val_Dof ;
	  }
	}
	else{		  
	  for (j = 0 ; j < Nbr_Dof ; j++){
	    for (k = 0 ; k < Current.NbrHar ; k+=2) {
	      if (Type_Quantity == QUANTITY_BF) {
		Val_Dof_r = (QuantityStorage_P->BasisFunction[j].Dof->Entity ==
			     Current.SubRegion)? 1. : 0. ;
		Val_Dof_i = 0. ;
	      }
	      else {
		Dof_GetComplexDofValue
		  (QuantityStorage_P->FunctionSpace->DofData,
		   QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT,
		   &Val_Dof_r, &Val_Dof_i) ;
	      }
	      Val[MAX_DIM*k  ]     += vBFxDof[j].Val[0] * Val_Dof_r ;
	      Val[MAX_DIM*k+1]     += vBFxDof[j].Val[1] * Val_Dof_r ;
	      Val[MAX_DIM*k+2]     += vBFxDof[j].Val[2] * Val_Dof_r ;
	      Val[MAX_DIM*(k+1)  ] += vBFxDof[j].Val[0] * Val_Dof_i ;
	      Val[MAX_DIM*(k+1)+1] += vBFxDof[j].Val[1] * Val_Dof_i ;
	      Val[MAX_DIM*(k+1)+2] += vBFxDof[j].Val[2] * Val_Dof_i ;
	    }
	  }
	}
	*Type_Value = VECTOR ;
	break ;
	
      default :
	Msg(GERROR, "Unknown Form type in 'Pos_FemInterpolation'");
	break;
      }
      
    }
    
    /* --------
       Integral
       -------- */

    /*
      Ceci, c'est nul a chier. Ce qu'il faut faire, c'est ne pas
      reinterpoler ici, mais laisser au Cal_Quantity dans
      Cal_IntegralQuantity le soin de reinterpoler directment la
      quantity local intervenant ds la qte integrale s'il y a lieu !

      Mais, comment faire avec l'integration analytique ?
    */

    
    else { 
      if (IQA.IntegrationCase_P->Type == ANALYTIC)	
	Cal_AnalyticIntegralQuantity (Current.Element = TheElement_P, 
				      QuantityStorage_P, Nbr_Dof, 
				      (void (**)())xFunctionBF, vBFxDof) ;
      else
	Cal_NumericalIntegralQuantity (Current.Element = TheElement_P, &IQA, 
				       QuantityStorage_P0, QuantityStorage_P, 
				       SubType_DefineQuantity, Nbr_Dof, 
				       (void (**)())xFunctionBF, vBFxDof) ;

      Type_Form = vBFxDof[0].Type ;
      
      /* interpolate (vBFxDof can be complex-valued) */
      
      if(SubType_DefineQuantity == NODOF){
	
	switch (Type_Form) {
	  
	case FORM0 :  case FORM3 :  case FORM3P :  case SCALAR :
	  for (k = 0 ; k < Current.NbrHar ; k++) 
	    Val[MAX_DIM*k] += vBFxDof[0].Val[MAX_DIM*k] ;
	  *Type_Value = SCALAR ;
	  break ;
	  
	case FORM1  :  case FORM1P :  case FORM2  :  case FORM2P :
	case FORM1S :  case FORM2S :
	case VECTOR :  case VECTORP :
	  for (k = 0 ; k < Current.NbrHar ; k++) {
	    Val[MAX_DIM*k]   += vBFxDof[0].Val[MAX_DIM*k] ;
	    Val[MAX_DIM*k+1] += vBFxDof[0].Val[MAX_DIM*k+1] ;
	    Val[MAX_DIM*k+2] += vBFxDof[0].Val[MAX_DIM*k+2] ;
	  }
	  *Type_Value = VECTOR ;
	  break ;
	  
	default :
	  Msg(GERROR, "Unknown Form type in 'Pos_FemInterpolation'");
	  break;
	}
	
      }
      else{
	
	switch (Type_Form) {
	  
	case FORM0 :  case FORM3 :  case FORM3P :  case SCALAR :
	  if(Current.NbrHar==1){
	    for (j = 0 ; j < Nbr_Dof ; j++){
	      Dof_GetRealDofValue
		(QuantityStorage_P->FunctionSpace->DofData,
		 QuantityStorage_P->BasisFunction[j].Dof,
		 &Val_Dof) ;
	      Val[0] += vBFxDof[j].Val[0] * Val_Dof ;
	    }
	  }
	  else{
	    for (j = 0 ; j < Nbr_Dof ; j++){
	      for (k = 0 ; k < Current.NbrHar ; k+=2) {
		Dof_GetComplexDofValue
		  (QuantityStorage_P->FunctionSpace->DofData,
		   QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT,
		   &Val_Dof_r, &Val_Dof_i) ;
		Val[MAX_DIM*k]     += 
		  vBFxDof[j].Val[MAX_DIM*k]     * Val_Dof_r -
		  vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_i ;
		Val[MAX_DIM*(k+1)] += 
		  vBFxDof[j].Val[MAX_DIM*k]     * Val_Dof_i +
		  vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_r ;
	      }
	    }
	  }
	  *Type_Value = SCALAR ;
	  break ;
	  
	case FORM1  :  case FORM1P :  case FORM2  :  case FORM2P :
	case FORM1S :  case FORM2S :
	case VECTOR :  case VECTORP :
	  if(Current.NbrHar==1){
	    for (j = 0 ; j < Nbr_Dof ; j++){
	      Dof_GetRealDofValue
		(QuantityStorage_P->FunctionSpace->DofData,
		 QuantityStorage_P->BasisFunction[j].Dof,
		 &Val_Dof) ;
	      Val[0] += vBFxDof[j].Val[0] * Val_Dof ;
	      Val[1] += vBFxDof[j].Val[1] * Val_Dof ;
	      Val[2] += vBFxDof[j].Val[2] * Val_Dof ;
	    }
	  }
	  else{	  
	    for (j = 0 ; j < Nbr_Dof ; j++){
	      for (k = 0 ; k < Current.NbrHar ; k+=2) {
		Dof_GetComplexDofValue
		  (QuantityStorage_P->FunctionSpace->DofData,
		   QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT,
		   &Val_Dof_r, &Val_Dof_i) ;
		Val[MAX_DIM*k]     += 
		  vBFxDof[j].Val[MAX_DIM*k]     * Val_Dof_r -
		  vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_i ;
		Val[MAX_DIM*(k+1)] += 
		  vBFxDof[j].Val[MAX_DIM*k]     * Val_Dof_i +
		  vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_r ;
		Val[MAX_DIM*k+1]     += 
		  vBFxDof[j].Val[MAX_DIM*k+1]     * Val_Dof_r -
		  vBFxDof[j].Val[MAX_DIM*(k+1)+1] * Val_Dof_i ;
		Val[MAX_DIM*(k+1)+1] += 
		  vBFxDof[j].Val[MAX_DIM*k+1]     * Val_Dof_i +
		  vBFxDof[j].Val[MAX_DIM*(k+1)+1] * Val_Dof_r ;
		Val[MAX_DIM*k+2]     += 
		  vBFxDof[j].Val[MAX_DIM*k+2]     * Val_Dof_r -
		  vBFxDof[j].Val[MAX_DIM*(k+1)+2] * Val_Dof_i ;
		Val[MAX_DIM*(k+1)+2] += 
		  vBFxDof[j].Val[MAX_DIM*k+2]     * Val_Dof_i +
		  vBFxDof[j].Val[MAX_DIM*(k+1)+2] * Val_Dof_r ;
	      }
	    }
	  }
	  *Type_Value = VECTOR ;
	  break ;
	  
	default :
	  Msg(GERROR, "Unknown Form type in 'Pos_FemInterpolation'");
	  break;
	}
	
      }
      
    }
    
    
    if (Type_DefineQuantity != INTEGRALQUANTITY)  break ;
    
  }  /* while (1) ... */
  
  if(UseNewGeo){
    GeoData_P = (struct GeoData *)List_Pointer(GeoData_L, GeoDataNum);
    Geo_SetCurrentGeoData(Current.GeoData = GeoData_P) ;
  }

  GetDP_End ;
}