fmm_caldtamatrices.c

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

  Msg(RESOURCES, "Before Disaggregation ");
  Msg(INFO, "Creation Disaggregation matrices");

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  NbrDir = Current.FMM.N ; 
  NbrHar = Current.NbrHar ;

  Current.FMM.Type = SCALAR ;

  Current.NbrHar = 2 ; /* FMM requires complex numbers, even with Laplace equation */

  for ( i_FMMEqu = 0 ; i_FMMEqu < NbrFMMEqu ; i_FMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + i_FMMEqu ;
    FMMObs = FMMmat_P->Obs ;
        
    EquationTerm_P = EquationTerm_P0 + FMMmat_P->EquTermIndex ;    
    FI = EquationTerm_P->Case.LocalTerm.Active ;

    if ( FMMmat_P->NbrCom == 0 )
      FMMmat_P->NbrCom = (Get_ValueFromForm(FI->Type_FormEqu)== VECTOR) ? 3 : 1 ;

    QuantityStorageEqu_P = FI->QuantityStorageEqu_P ;
    QuantityStorageDof_P = FI->QuantityStorageDof_P ;

    if( EquationTerm_P->Case.LocalTerm.FMMIntegrationMethodIndex == -1 )
      EquationTerm_P->Case.LocalTerm.FMMIntegrationMethodIndex = 
	EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex ;

    FMMIntegrationCaseEqu_L = 
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod,
		    EquationTerm_P->Case.LocalTerm.FMMIntegrationMethodIndex))
      ->IntegrationCase ;

    FMMCriterionIndexEqu = 
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod,
		    EquationTerm_P->Case.LocalTerm.FMMIntegrationMethodIndex))
      ->CriterionIndex ;

    List_Read(Problem_S.FMMGroup, FMMObs, &FMMGroup_S ) ;   
    Nbr_Group  = List_Nbr(FMMGroup_S.List ) ;
   
    FMMData_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;   

    for ( i_Group = 0 ; i_Group < Nbr_Group ; i_Group++ ) {
      FMMData_P =  FMMData_P0 + i_Group ;
      Nbr_ElmsInGroup = List_Nbr(FMMData_P->Element);

      ElmtsGr = (int*)(FMMData_P->Element->array) ;

      List_Read(FMMmat_P->D_L, i_Group, &Disagg_M) ;
      List_Read(FMMmat_P->NumEqu, i_Group, &NumEqu_L) ;
      Nbr_EquList = List_Nbr(NumEqu_L) ;

      /* Initializing Disaggregation matrices */
      if (Disagg_M == NULL){
	Disagg_M = (double**)Malloc(Nbr_EquList*sizeof(double*)) ;
	for (j = 0 ; j < Nbr_EquList ; j++)
	  Disagg_M[j] = (double*)Calloc(2 * FMMmat_P->NbrCom * NbrDir, sizeof(double)) ;
      }

      Current.FMM.Xgc  = FMMData_P->Xgc ; /* FMM group center */
      Current.FMM.Ygc  = FMMData_P->Ygc ;
      Current.FMM.Zgc  = FMMData_P->Zgc ;
      Current.FMM.Robs = FMMData_P->Rmax ;

      for (i_Element = 0 ; i_Element < Nbr_ElmsInGroup ; i_Element++) {
	Element.GeoElement = Geo_GetGeoElementOfNum(ElmtsGr[i_Element]) ;
	Element.Num    = Element.GeoElement->Num ;
	Element.Type   = Element.GeoElement->Type ;
	Current.Region = Element.Region = Element.GeoElement->Region ;
	if (Element.GeoElement->FMMGroup != i_Group) 
	  Msg(GERROR, "Element.GeoElement->FMMGroup != i_Group");
	Element.FMMGroup = i_Group ;
	
	QuantityStorageEqu_P->NumLastElementForFunctionSpace = Element.Num ;
	Get_DofOfElement
	  (&Element, QuantityStorageEqu_P->FunctionSpace, QuantityStorageEqu_P,
 	   NULL) ;

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

	Current.Element = &Element ; 
	Element.ElementSource = &Element ;
	Current.ElementSource = Element.ElementSource ;

	FI->xChangeOfCoordinatesEqu = 
	  (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormEqu) ;

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

	IntegrationCase_P =
	  Get_IntegrationCase(&Element, FMMIntegrationCaseEqu_L, FMMCriterionIndexEqu);

	switch (IntegrationCase_P->Type) {
	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.FMMIntegrationMethodIndex))->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) ;
	    }/* ChangeCoord */
	    
	    Current.x = Current.y = Current.z = 0. ;

	    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] ;
	    }
	
	    Factor = (FI->Flag_ChangeCoord) ? weight * fabs(Element.DetJac) : weight ;

	    Current.FMM.Flag_GF = FMM_DISAGGREGATION ;
	    ((CASTF2V)QuantityStorageDof_P->DefineQuantity->IntegralQuantity.FunctionForFMM.Fct)
	      (&QuantityStorageDof_P->DefineQuantity->IntegralQuantity.FunctionForFMM,
	       GFValue, GFValue ) ;
	    
	    Cal_ZeroValue(&TermFct);
	    Cal_TermsforDisaggregation( EquationTerm_P, QuantityStorage_P0, &TermFct ) ;

	    /* Test Functions */		

	    for (i = 0 ; i < Nbr_Equ ; i++) {
	      if (QuantityStorageEqu_P->BasisFunction[i].Dof->Type == DOF_UNKNOWN){
		
		xFunctionBFEqu[i]
		  (&Element,
		   QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement+1,
		   Current.u, Current.v, Current.w, vBFuEqu) ;
		vBFxEqu[0] = vBFxEqu[1]= vBFxEqu[2] = 0. ;/* Initialization */		
		((void (*)(struct Element*, double*, double*))
		 FI->xChangeOfCoordinatesEqu) (&Element, vBFuEqu, vBFxEqu) ;
		
		vBFxEquV.Type = Get_ValueFromForm(FI->Type_FormEqu) ;
		vBFxEquV.Val[0]         = vBFxEqu[0] * Factor ;
		vBFxEquV.Val[1]         = vBFxEqu[1] * Factor ;
		vBFxEquV.Val[2]         = vBFxEqu[2] * Factor ;
		vBFxEquV.Val[MAX_DIM]   = 0. ;
		vBFxEquV.Val[MAX_DIM+1] = 0. ;
		vBFxEquV.Val[MAX_DIM+2] = 0. ;
		
		switch(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity ){
		case CWQ_EXP_TIME_DOF :
		case CWQ_FCT_TIME_DOF :
		  Cal_ProductValue(&TermFct, &vBFxEquV, &vBFxEquV) ;
		  break ;
		case CWQ_FCT_PVEC_DOF :
		  Cal_CrossProductValue(&TermFct, &vBFxEquV, &vBFxEquV) ;
		  break ;
		default:
		  break ;
		}
		
		for (i_FMM = 0 ; i_FMM < NbrDir ; i_FMM++) 		
		  Cal_ProductValue(&vBFxEquV, &GFValue[i_FMM], &BFGFValue[i_FMM]);
	   
		j = List_ISearch(NumEqu_L, 
				 &QuantityStorageEqu_P->BasisFunction[i].Dof->
				 Case.Unknown.NumDof, fcmp_int) ;
		
		if (j != -1) Cal_AddValueArray2DoubleArray(BFGFValue, Disagg_M[j], NbrDir) ; 
		else  Msg(GERROR, "Wrong NumEqu %d for Disagg",
			  QuantityStorageEqu_P->BasisFunction[i].Dof->Case.Unknown.NumDof) ;
		
	      }/* if DOF_UNKNOWN */      
	    } /* for i Nbr_Equ */
	  } /* for i_IntPoint ... */            
	  break ; /* case GAUSS */

	case ANALYTIC :
	  Current.FMM.Flag_GF = FMM_DISAGGREGATION ;
	  Cal_ZeroValue(&TermFct);
	  Cal_TermsforDisaggregation( EquationTerm_P, QuantityStorage_P0, &TermFct ) ;
	  GF_FMMLaplacexForm( &Element, QuantityStorageEqu_P, Nbr_Equ, (void (*)())xFunctionBFEqu, NumEqu_L, TermFct, Disagg_M );
	  break;
	}/*switch Integration_Case */	
      } /* i_Element */
      List_Write(FMMmat_P->D_L, i_Group, &Disagg_M) ;
    }/* i_Group */
 
  }/* i_FMMEqu */

  Current.FMM.Flag_GF = FMM_DIRECT ;
  Current.NbrHar = NbrHar ;

  GetDP_End ;
}


/*-------------------------------------------------------------------------------*/
/*           C a l _ T e r m s f o r D i s a g g r e g a t i o n                 */
/*-------------------------------------------------------------------------------*/

void Cal_TermsforDisaggregation( struct EquationTerm * EquationTerm_P,
				 struct QuantityStorage * QuantityStorage_P0,
				 struct Value *TermFct ){

  GetDP_Begin("Cal_TermsforDisaggregation") ;

  switch(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity){
  case CWQ_EXP_TIME_DOF:
    Get_ValueOfExpression
      (Problem_Expression0 +
       EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical,
       QuantityStorage_P0, 0., 0., 0., 
       TermFct) ;
    break;
  case CWQ_FCT_TIME_DOF:
  case CWQ_FCT_PVEC_DOF:
    ((CASTF2V)EquationTerm_P->Case.LocalTerm.Term.FunctionForCanonical.Fct)
      (&EquationTerm_P->Case.LocalTerm.Term.FunctionForCanonical, TermFct, TermFct) ;
    break;
  default: break ; 
  }

  GetDP_End ;

}


/*-------------------------------------------------------------------------------*/
/*           C a l _ F M M G a l e r k i n A g g r e g a t i o n           */
/*-------------------------------------------------------------------------------*/

void  Cal_FMMGalerkinAggregation(struct EquationTerm     * EquationTerm_P0,
				 struct QuantityStorage  * QuantityStorage_P0) {
  
  struct EquationTerm       * EquationTerm_P ;
  struct FemLocalTermActive * FI ;
  struct QuantityStorage    * QuantityStorageDof_P ;
  struct IntegrationCase    * IntegrationCase_P ;
  struct Quadrature         * Quadrature_P ;
  struct Value                vBFxDofV ;
  struct Value                Val0 ;
  struct Value                GFValue [NBR_MAX_DIR],  BFGFValue [NBR_MAX_DIR] ;
  struct Element              Element ;
  struct FMMData            * FMMData_P0, * FMMData_P ; 
  struct FMMGroup             FMMGroup_S ;  
  struct FMMmat             * FMMmat_P0 , * FMMmat_P ;
  double                    ** Aggreg_M ;
  List_T                    * NumDof_L ;

  int     Nbr_Dof ;
  int     Nbr_Group, i_Group, Nbr_ElmsInGroup, i_Element, *ElmtsGr ;
  int     FMMSrc, NbrFMMEqu, i_FMMEqu, NbrDir, i_FMM ;
  int     i, j, Type_Dimension, Nbr_IntPoints, i_IntPoint ;
  int     Nbr_DofList, NbrHar ;
  
  double  weight, Factor = 1. ;
  double  vBFuDof [MAX_DIM], vBFxDof [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*);
  void   (*FunctionProd)() ;

  List_T * FMMIntegrationCaseDof_L ;