cal_galerkintermoffemequation.c

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

  /*  ----------------------------------------------------------------------  */

  if (!(Nbr_Equ = QuantityStorageEqu_P->NbrElementaryBasisFunction)) {
    GetDP_End ;
  }

  Get_FunctionValue(Nbr_Equ, (void (**)())xFunctionBFEqu,
		    EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu,
		    QuantityStorageEqu_P, &FI->Type_FormEqu) ;

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

  if (FI->SymmetricalMatrix){
    Nbr_Dof = Nbr_Equ ;
  }
  else{
    switch (FI->Type_DefineQuantityDof) {
    case LOCALQUANTITY :
      Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ;
      Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
			EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
			QuantityStorageDof_P, &FI->Type_FormDof) ;
      break ;
    case INTEGRALQUANTITY :
      Get_InitElementSource(Element,
			    QuantityStorageDof_P->DefineQuantity->IntegralQuantity.InIndex) ;  
      break ;
    case NODOF : 
      Nbr_Dof = 1 ;  
      break ;
    }
  }

  /*  -------------------------------------------------------------------  */
  /*  G e t   I n t e g r a t i o n   M e t h o d                          */
  /*  -------------------------------------------------------------------  */

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

  /*  -------------------------------------------------------------------  */
  /*  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) ;

  /*  ------------------------------------------------------------------------  */
  /*  ------------------------------------------------------------------------  */
  /*  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           */
  /*  ------------------------------------------------------------------------  */
  /*  ------------------------------------------------------------------------  */


  /* Loop on source elements (> 1 only if integral quantity) */
  
  while (1) {  

    if (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) {

      if (!Flag_FMM)
	NextElement = Get_NextElementSource(Element->ElementSource) ;
      else{
	NextElement = Get_NextElementSourceNeighbour(Element) ;
      }
      
      if (NextElement) {
	/* Get DOF of source element */
      
	Get_DofOfElement(Element->ElementSource,
			 QuantityStorageDof_P->FunctionSpace,
			 QuantityStorageDof_P, NULL) ;
	
	/* Get function value for shape function */
      
	Get_NodesCoordinatesOfElement(Element->ElementSource) ;
	Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ;
	Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
			  QuantityStorageDof_P->DefineQuantity->IntegralQuantity.TypeOperatorDof, 
			  QuantityStorageDof_P, &FI->IntegralQuantityActive.Type_FormDof) ;
	
	/* Initialize Integral Quantity (integration & jacobian) */
	
	Cal_InitIntegralQuantity(Element, &FI->IntegralQuantityActive, 
				   QuantityStorageDof_P);
      }
      else { 
	break ;	  
      } /* if NextElement */
    } /* if INTEGRALQUANTITY */
      
   
    if (FI->SymmetricalMatrix)
      for (i = 0 ; i < Nbr_Equ ; i++)  for (j = 0 ; j <= i      ; j++) 
	for (k = 0 ; k < Current.NbrHar ; k++)  Ek[i][j][k] = 0. ;
    else
      for (i = 0 ; i < Nbr_Equ ; i++)  for (j = 0 ; j < Nbr_Dof ; j++)
	for (k = 0 ; k < Current.NbrHar ; k++)  Ek[i][j][k] = 0. ;
    
   
    switch (IntegrationCase_P->Type) {
      
      /*  -------------------------------------  */
      /*  Q U A D R A T U R E                    */
      /*  -------------------------------------  */
      
    case GAUSS :  
    case GAUSSLEGENDRE :
     
      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 ;
           
      
      for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
	
	Get_IntPoint(Nbr_IntPoints, i_IntPoint, 
		     &Current.u, &Current.v, &Current.w, &weight) ;
	
	if (FI->Flag_ChangeCoord) {
	  Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ;

	  Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
	  
	  if (FI->Flag_InvJac)
	    Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
			      &Element->Jac, &Element->InvJac) ;
	}
	
	/* Test Functions */
	
	if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE)
	  Get_ValueOfExpressionByIndex
	    (EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical_Equ,
	     QuantityStorage_P0, Current.u, Current.v, Current.w, &CanonicExpression_Equ) ;
	
	for (i = 0 ; i < Nbr_Equ ; i++) {
	  xFunctionBFEqu[i]
	    (Element,
	     QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement+1,
	     Current.u, Current.v, Current.w, vBFuEqu[i]) ;
	  ((void (*)(struct Element*, double*, double*))
	   FI->xChangeOfCoordinatesEqu) (Element, vBFuEqu[i], vBFxEqu[i]) ;
	  
	  
	  if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE){
	    V1.Type = Get_ValueFromForm(FI->Type_FormEqu);
	    V1.Val[0]         = vBFxEqu[i][0] ;
	    V1.Val[1]         = vBFxEqu[i][1] ;
	    V1.Val[2]         = vBFxEqu[i][2] ;
	    V1.Val[MAX_DIM]   = 0;
	    V1.Val[MAX_DIM+1] = 0;
	    V1.Val[MAX_DIM+2] = 0;
	    
	    switch(EquationTerm_P->Case.LocalTerm.Term.OperatorTypeForCanonical_Equ){
	    case OP_TIME : Cal_ProductValue (&CanonicExpression_Equ,&V1,&V2); break;
	    case OP_CROSSPRODUCT : Cal_CrossProductValue (&CanonicExpression_Equ,&V1,&V2); break;
	    default : Msg(GERROR, "Unknown operation in Equation");
	    }
	    vBFxEqu[i][0] = V2.Val[0];
	    vBFxEqu[i][1] = V2.Val[1];
	    vBFxEqu[i][2] = V2.Val[2];
	  }
	  
	} /* for Nbr_Equ */
	
	
	/* Shape Functions (+ surrounding expression) */
	
	Current.Element = Element ;
	Cal_vBFxDof(EquationTerm_P, FI, 
		    QuantityStorage_P0, QuantityStorageDof_P,
		    Nbr_Dof, xFunctionBFDof, vBFxEqu, vBFxDof);

	Factor = (FI->Flag_ChangeCoord) ? weight * fabs(Element->DetJac) : weight ;

	/* Product and assembly in elementary submatrix             (k?-1.:1.)*   */

	if (FI->SymmetricalMatrix)
	  for (i = 0 ; i < Nbr_Equ ; i++)  for (j = 0 ; j <= i ; j++)
	    for (k = 0 ; k < Current.NbrHar ; k++)
	      Ek[i][j][k] += Factor *
		((double (*)(double*, double*))
		 FI->Cal_Productx) (vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM*k])) ;
	else
	  for (i = 0 ; i < Nbr_Equ ; i++)  for (j = 0 ; j < Nbr_Dof ; j++)
	    for (k = 0 ; k < Current.NbrHar ; k++)
	      Ek[i][j][k] += Factor *
		((double (*)(double*, double*))
		 FI->Cal_Productx) (vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM*k])) ;

      }  /* for i_IntPoint ... */            
      break ; /* case GAUSS */
      
      
      /*  -------------------------------------  */
      /*  A N A L Y T I C                        */
      /*  -------------------------------------  */
      
    case ANALYTIC :
      
      if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity ==
	  CWQ_DOF) {
	Factor = 1. ;
      }
      if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity ==
	  CWQ_EXP_TIME_DOF) {
	if (EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical >= 0) {
	  Get_ValueOfExpression
	    (Problem_Expression0 +
	     EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical,
	     QuantityStorage_P0, 0., 0., 0., &CoefPhys) ;
	  Factor = CoefPhys.Val[0] ;
	}
      }
      else {
	Msg(GERROR, "Bad expression for full analytic integration");
      }
      
      if (FI->SymmetricalMatrix) {
	for (i = 0 ; i < Nbr_Equ ; i++)  for (j = 0 ; j <= i ; j++)
	  Ek[i][j][0] = Factor *
	    Cal_AnalyticIntegration
	    (Element, (void (*)())xFunctionBFEqu[i], (void (*)())xFunctionBFEqu[j], i, j,
	     FI->Cal_Productx) ;
      }
      else {
	switch (FI->Type_DefineQuantityDof) {
	case LOCALQUANTITY :
	  for (i = 0 ; i < Nbr_Equ ; i++)  for (j = 0 ; j < Nbr_Dof ; j++)
	    Ek[i][j][0] = Factor *
	      Cal_AnalyticIntegration
	      (Element, (void (*)())xFunctionBFEqu[i], (void (*)())xFunctionBFDof[j], i, j,
	       FI->Cal_Productx) ;
	  break;
	default :
	  Msg(GERROR, "Exterior analytical integration not implemented");
	  break;
	}
      }
      
      break ; /* case ANALYTIC */
      
    default :
      Msg(GERROR, "Unknown type of Integration method (%s)",
	  ((struct IntegrationMethod *)
	   List_Pointer(Problem_S.IntegrationMethod,
			EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
      break;
      
    }

    /* Complete elementary matrix if symmetrical */
    /* ----------------------------------------- */
    
    if (FI->SymmetricalMatrix)
      for (i = 1 ; i < Nbr_Equ ; i++)  
	for (j = 0 ; j < i ; j++)
	  for (k = 0 ; k < Current.NbrHar ; k++)  
	    Ek[j][i][k] = Ek[i][j][k] ;
    
    if(Flag_VERBOSE == 10) {
      printf("Galerkin = ") ;
      for (j = 0 ; j < Nbr_Dof ; j++)
	Print_DofNumber(QuantityStorageDof_P->BasisFunction[j].Dof) ; 
      printf("\n") ;
      for (i = 0 ; i < Nbr_Equ ; i++) { 
	Print_DofNumber(QuantityStorageEqu_P->BasisFunction[i].Dof) ; 
	printf("[ ") ;
	for (j = 0 ; j < Nbr_Dof ; j++) {
	  printf("(") ;
	  for(k = 0 ; k < Current.NbrHar ; k++) printf("% .5e ", Ek[i][j][k]) ;
	  printf(")") ;
	}
	printf("]\n") ;
      }
    }
    
    /* Assembly in global matrix */
    /* ------------------------- */
   
    for (i = 0 ; i < Nbr_Equ ; i++)  
      for (j = 0 ; j < Nbr_Dof ; j++)
	((void (*)(struct Dof*, struct Dof*, double*)) 
	 FI->Function_AssembleTerm)
	  (QuantityStorageEqu_P->BasisFunction[i].Dof,
	   QuantityStorageDof_P->BasisFunction[j].Dof,  Ek[i][j]) ;
    
    /* Exit while(1) directly if not integral quantity */

    if (FI->Type_DefineQuantityDof != INTEGRALQUANTITY)  break ;

  }  /* while (1) ... */
   
  GetDP_End ;
}