f_multihar.c

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

  struct FemLocalTermActive  * FI ;
  List_T * WholeQuantity_L;
  struct WholeQuantity   *WholeQuantity_P0 ;
  int i_WQ ;

  GetDP_Begin("Cal_InitGalerkinTermOfFemEquation_MHJacNL");

  FI = EquationTerm_P->Case.LocalTerm.Active ;
  FI->MHJacNL = 0 ;

  /* search for MHJacNL-term(s) */ 
  WholeQuantity_L = EquationTerm_P->Case.LocalTerm.Term.WholeQuantity ;
  WholeQuantity_P0 = (struct WholeQuantity*)List_Pointer(WholeQuantity_L, 0) ;
  i_WQ = 0 ; while ( i_WQ < List_Nbr(WholeQuantity_L) &&
		     (WholeQuantity_P0 + i_WQ)->Type != WQ_MHJACNL) i_WQ++ ;

  if (i_WQ == List_Nbr(WholeQuantity_L) ) 
    GetDP_End ;   /* no MHJacNL stuff, let's get the hell out of here ! */


  /* check if Galerkin term produces symmetrical contribution to system matrix */
  if (!FI->SymmetricalMatrix) 
     Msg(GERROR, "Galerkin term with MHJacNL must be symmetrical");
  
  if (List_Nbr(WholeQuantity_L) == 3){
    if (i_WQ != 0 || 
	EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity != 1 ||
	(WholeQuantity_P0 + 2)->Type != WQ_BINARYOPERATOR ||
	(WholeQuantity_P0 + 2)->Case.Operator.TypeOperator != OP_TIME)
      Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL");
    FI->MHJacNL_Factor = 1.; 
  } else if (List_Nbr(WholeQuantity_L) == 5){
    if ((WholeQuantity_P0 + 0)->Type != WQ_CONSTANT ||
	i_WQ != 1 || 
	(WholeQuantity_P0 + 2)->Type != WQ_BINARYOPERATOR ||
	(WholeQuantity_P0 + 2)->Case.Operator.TypeOperator != OP_TIME ||
	EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity != 3 ||
	(WholeQuantity_P0 + 4)->Type != WQ_BINARYOPERATOR ||
	(WholeQuantity_P0 + 4)->Case.Operator.TypeOperator != OP_TIME)
      Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL");
    FI->MHJacNL_Factor = WholeQuantity_P0->Case.Constant ;
    /* printf(" Factor = %e \n" , FI->MHJacNL_Factor); */
  } else {
    Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL (%d terms) ", 
	List_Nbr(WholeQuantity_L));
  }

  if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE) 
    Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL");
    
  if (EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative != JACNL_)
    Msg(GERROR, "MHJacNL can only be used with JACNL") ;

  FI->MHJacNL = 1 ;
  FI->MHJacNL_Index  = (WholeQuantity_P0 + i_WQ)->Case.MHJacNL.Index ;

  FI->MHJacNL_HarOffSet = 2 * (WholeQuantity_P0 + i_WQ)->Case.MHJacNL.FreqOffSet ;  
  if (FI->MHJacNL_HarOffSet > Current.NbrHar-2){ 
    Msg(WARNING, "FreqOffSet in 'MHJacNL' cannot exceed %d => set to %d ",
	Current.NbrHar/2-1, Current.NbrHar/2-1) ;
    FI->MHJacNL_HarOffSet = Current.NbrHar-2 ;
  }

  MH_Get_InitData(2, (WholeQuantity_P0 + i_WQ)->Case.MHJacNL.NbrPoints,
		  &FI->MHJacNL_NbrPointsX, &FI->MHJacNL_H, &FI->MHJacNL_HH,
		  &FI->MHJacNL_t, &FI->MHJacNL_w) ;

  GetDP_End ;
}


#define OFFSET (iHar < NbrHar-OffSet)? 0 : iHar-NbrHar+OffSet+2-iHar%2


/* --------------------------------------------------------------------------- */
/*  C a l _ G a l e r k i n T e r m O f F e m E q u a t i o n _ M H J a c N L  */
/* --------------------------------------------------------------------------- */

void  Cal_GalerkinTermOfFemEquation_MHJacNL(struct Element          * Element,
					    struct EquationTerm     * EquationTerm_P,
					    struct QuantityStorage  * QuantityStorage_P0) {
  
  struct FemLocalTermActive * FI ;
  struct QuantityStorage    * QuantityStorage_P;
  struct IntegrationCase    * IntegrationCase_P ;
  struct Quadrature         * Quadrature_P ;

  int     Nbr_Dof, NbrHar ;
  double  vBFuDof [NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
  double  vBFxDof [NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
  double  Val_Dof [NBR_MAX_BASISFUNCTIONS][NBR_MAX_HARMONIC] ;
  double  Val     [3*NBR_MAX_BASISFUNCTIONS];

  int     i, j, k, Type_Dimension,  Nbr_IntPoints, i_IntPoint ;
  int     iTime, iDof, jDof, iHar, jHar, nVal1, nVal2 = 0, iVal1, iVal2, Type1;
  double **H, ***HH, Factor, plus, plus0, weightIntPoint;
  int NbrPointsX, OffSet;
  struct Expression * Expression_P;
  struct Dof * Dofi, *Dofj;
  struct Value t_Value;
  gMatrix * Jac;

  double one=1.0 ;
  int iPul, ZeroHarmonic, DcHarmonic;

  double E_MH[NBR_MAX_BASISFUNCTIONS][NBR_MAX_BASISFUNCTIONS][NBR_MAX_HARMONIC][NBR_MAX_HARMONIC];
  double E_D[NBR_MAX_HARMONIC][NBR_MAX_HARMONIC][MAX_DIM];

  void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
    (struct Element * Element, int NumEntity, 
     double u, double v, double w, double Value[] ) ;
  double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
  void (*Get_IntPoint)(int,int,double*,double*,double*,double*);
  double (*Get_Product)(double*,double*,double*) = 0;

  /* static double eps; */
 
  GetDP_Begin("Cal_GalerkinTermOfFemEquation_MHJacNL");

  FI = EquationTerm_P->Case.LocalTerm.Active ;
  QuantityStorage_P = FI->QuantityStorageDof_P ;

  /*  ----------------------------------------------------------------------  */
  /*  G e t   F u n c t i o n V a l u e  f o r  t e s t  f u n c t i o n s    */
  /*  ----------------------------------------------------------------------  */

  if (!(Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction)) GetDP_End ;

  Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
		    EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
		    QuantityStorage_P, &FI->Type_FormDof) ;

  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[j][k], &Val_Dof[j][k+1]) ;
  
  /* printf("Type1 = %d\n",  FI->Type_FormDof); */
  nVal1 = NbrValues_Type (Type1 = Get_ValueFromForm(FI->Type_FormDof)) ;



  /*  -------------------------------------------------------------------  */
  /*  G e t   J a c o b i a n   M e t h o d                                */
  /*  -------------------------------------------------------------------  */
  
  i = 0 ; while ((i < FI->NbrJacobianCase) &&
		 ((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) &&
		 !List_Search
		 (((struct Group *)List_Pointer(Problem_S.Group, j))
		  ->InitialList, &Element->Region, fcmp_int) ) i++ ;

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

  Element->JacobianCase = FI->JacobianCase_P0 + i ;

  Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
    Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
			 Element->Type, &Type_Dimension) ;

  if (FI->Flag_ChangeCoord) Get_NodesCoordinatesOfElement(Element) ;

 

  /*  integration in space  */

  IntegrationCase_P =
    Get_IntegrationCase(Element, FI->IntegrationCase_L, FI->CriterionIndex);

  switch (IntegrationCase_P->Type) {    
  case ANALYTIC :
    Msg(GERROR, "Analytical integration not implemented for MHJacNL");
  }  

  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)",
	Get_StringForDefine(Element_Type, Element->Type),
	((struct IntegrationMethod *)
	 List_Pointer(Problem_S.IntegrationMethod,
			EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
      
  Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
  Get_IntPoint  = (void(*)(int,int,double*,double*,double*,double*))
    Quadrature_P->Function ;


  /*  integration in fundamental time period  */

  NbrPointsX   = FI->MHJacNL_NbrPointsX;
  HH           = FI->MHJacNL_HH;
  H            = FI->MHJacNL_H ;
  Expression_P = Problem_Expression0 + FI->MHJacNL_Index;
  OffSet       = FI->MHJacNL_HarOffSet;
  Factor       = FI->MHJacNL_Factor;



  /*  ------------------------------------------------------------------------  */
  /*  ------------------------------------------------------------------------  */
  /*  C o m p u t a t i o n   o f   E l e m e n t a r y   m a t r i x           */
  /*  ------------------------------------------------------------------------  */
  /*  ------------------------------------------------------------------------  */

  NbrHar = Current.NbrHar ;

  /* special treatment of DC-term and associated dummy sinus-term */
  DcHarmonic = NbrHar;
  ZeroHarmonic = 0;
  for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) 
    if (!Current.DofData->Val_Pulsation[iPul]){
      DcHarmonic = 2*iPul ;
      ZeroHarmonic = 2*iPul+1 ;
      break;
    }

  /* resetting elementary matrix */
  for (iDof = 0 ; iDof < Nbr_Dof ; iDof++)
    for (jDof = 0 ; jDof <= iDof ; jDof++)
      for (iHar = 0 ; iHar < NbrHar ; iHar++)
	for (jHar = OFFSET ; jHar <= iHar ; jHar++)
	  E_MH[iDof][jDof][iHar][jHar] = 0. ;
  
  /* volume integration over element */ 
  for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
    
    Get_IntPoint(Nbr_IntPoints, i_IntPoint, 
		 &Current.u, &Current.v, &Current.w, &weightIntPoint) ;
    
    if (FI->Flag_ChangeCoord) {
      Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ;
      
      Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
      weightIntPoint *= fabs(Element->DetJac) ;
      
      if (FI->Flag_InvJac)
	Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
			  &Element->Jac, &Element->InvJac) ;
    }
	
    /* Test and shape Functions (are the same) */
	
    for (i = 0 ; i < Nbr_Dof ; i++) {
      xFunctionBFDof[i]
	(Element,
	 QuantityStorage_P->BasisFunction[i].NumEntityInElement+1,
	 Current.u, Current.v, Current.w, vBFuDof[i]) ;
      ((void (*)(struct Element*, double*, double*))
       FI->xChangeOfCoordinatesEqu) (Element, vBFuDof[i], vBFxDof[i]) ;
    } 
	
    switch (Type1) {        
    case SCALAR :
      for (k = 0 ; k < NbrHar ; k++){ 
	Val[k] = 0.;
	for (j = 0 ; j < Nbr_Dof ; j++)
	  Val[k] += vBFxDof[j][0] * Val_Dof[j][k] ;
      }
      break ;
    case VECTOR :
      for (k = 0 ; k < NbrHar ; k++){
	Val[3*k] = Val[3*k+1] = Val[3*k+2] = 0.;
	for (j = 0 ; j < Nbr_Dof ; j++){
	  Val[3*k  ] += vBFxDof[j][0] * Val_Dof[j][k] ;
	  Val[3*k+1] += vBFxDof[j][1] * Val_Dof[j][k] ;
	  Val[3*k+2] += vBFxDof[j][2] * Val_Dof[j][k] ;
	}
      }
      break ;
    }

   
    Current.NbrHar = 1;  /* evaluation in time domain */

    /* time integration over fundamental period */
    for (iTime = 0 ; iTime < NbrPointsX ; iTime++) {  

      t_Value.Type = Type1 ;
      for (iVal1 = 0 ; iVal1 < nVal1 ; iVal1++){
	t_Value.Val[iVal1] = 0; 
	for (iHar = 0 ; iHar < NbrHar ; iHar++)
	  t_Value.Val[iVal1] += H[iTime][iHar] * Val[iHar*nVal1+iVal1] ;
      }      
      
      Get_ValueOfExpression(Expression_P, QuantityStorage_P0, 
			    Current.u, Current.v, Current.w, &t_Value);
      
      if (!iTime){
	if (!i_IntPoint){
	  nVal2 = NbrValues_Type (t_Value.Type) ; 
	  Get_Product = (double(*)(double*,double*,double*))
	    Get_RealProductFunction_Type1xType2xType1 (Type1, t_Value.Type) ;
	}
	for (iHar = 0 ; iHar < NbrHar ; iHar++)
	  for (jHar = OFFSET ; jHar <= iHar ; jHar++)
	    for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
	      E_D[iHar][jHar][iVal2] = 0. ;
      }

      for (iHar = 0 ; iHar < NbrHar ; iHar++)
	for (jHar = OFFSET  ; jHar <= iHar ; jHar++){
	  for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
	    E_D[iHar][jHar][iVal2] += HH[iTime][iHar][jHar] * t_Value.Val[iVal2] ;
	}

    } /* for iTime ... */

    /*
    if (!eps) {
      printf("enter value for eps\n");
      scanf("%lf",&eps);
      printf("eps = %f\n",eps);
    }


    for (iHar = 0 ; iHar < NbrHar ; iHar++)
      for (jHar = OFFSET  ; jHar <= iHar ; jHar++){
	for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
	  if ( E_D[iHar][jHar][iVal2] * E_D[iHar][jHar][iVal2] < 
	      eps * eps * fabs(E_D[iHar][iHar][iVal2] * E_D[jHar][jHar][iVal2]) )
	    E_D[iHar][jHar][iVal2]=0 ;
      }

    */

    for (iDof = 0 ; iDof < Nbr_Dof ; iDof++)
      for (jDof = 0 ; jDof <= iDof ; jDof++)
	for (iHar = 0 ; iHar < NbrHar ; iHar++)
	  for (jHar = OFFSET  ; jHar <= iHar ; jHar++){
	    E_MH[iDof][jDof][iHar][jHar] += weightIntPoint * 
	      Get_Product(vBFxDof[iDof], E_D[iHar][jHar], vBFxDof[jDof]) ;

	    /*  printf("%d %d %d %d  %e\n", iDof, jDof, iHar, jHar,
	        E_MH[iDof][jDof][iHar][jHar]) ; */
	  }

    Current.NbrHar = NbrHar ;
 
  }  /* for i_IntPoint ... */            


  /*  --------------------------------------------------------------------  */
  /*  A d d   c o n t r i b u t i o n   t o  J a c o b i a n   M a t r i x  */
  /*  --------------------------------------------------------------------  */


  Jac = &Current.DofData->Jac;

  for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
    Dofi = QuantityStorage_P->BasisFunction[iDof].Dof ;    
    for (jDof = 0 ; jDof <= iDof ; jDof++){
      Dofj = QuantityStorage_P->BasisFunction[jDof].Dof ;

      for (iHar = 0 ; iHar < NbrHar ; iHar++)
	for (jHar = OFFSET ; jHar <= iHar ; jHar++){
	  plus = plus0 = Factor * E_MH [iDof][jDof] [iHar][jHar] ;
	  if(jHar==DcHarmonic && iHar!=DcHarmonic) { plus0 *= 1. ; plus *= 2. ;}
	  Dof_AssembleInMat(Dofi+iHar, Dofj+jHar, 1, &plus, Jac, NULL) ;
	  if(iHar != jHar)
	      Dof_AssembleInMat(Dofi+jHar, Dofj+iHar, 1, &plus0, Jac, NULL) ;
	  if(iDof != jDof){
	      Dof_AssembleInMat(Dofj+iHar, Dofi+jHar, 1, &plus, Jac, NULL) ;
	      if(iHar != jHar)
		Dof_AssembleInMat(Dofj+jHar, Dofi+iHar, 1, &plus0, Jac, NULL) ;
	  }
	}
    }
  }


  /* dummy 1's on the diagonal for sinus-term of dc-component */

  if (ZeroHarmonic) {
    for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
      Dofi = QuantityStorage_P->BasisFunction[iDof].Dof + ZeroHarmonic ;
      Dof_AssembleInMat(Dofi, Dofi, 1, &one, Jac, NULL) ;
    }
  }
 

  GetDP_End ;
}


#undef OFFSET