fmm_caldtamatrices.c

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

  int      FMMCriterionIndexDof ;


  GetDP_Begin("Cal_FMMGalerkinAggregation");

  Msg(RESOURCES, "Before Aggregation ");
  Msg(INFO, "Creation Aggregation 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.NbrHar = 2 ; /* ALWAYS complex numbers in FMM data matrices*/

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

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

    QuantityStorageDof_P = FI->QuantityStorageDof_P ;

    if( FI->QuantityStorageDof_P->DefineQuantity->
	IntegralQuantity.FMMIntegrationMethodIndex == -1 )
      FI->QuantityStorageDof_P->DefineQuantity->IntegralQuantity.FMMIntegrationMethodIndex =
	FI->QuantityStorageDof_P->DefineQuantity->IntegralQuantity.IntegrationMethodIndex ;

    FMMIntegrationCaseDof_L= 
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod,
		    FI->QuantityStorageDof_P->DefineQuantity->
		    IntegralQuantity.FMMIntegrationMethodIndex)) ->IntegrationCase ;
      
    FMMCriterionIndexDof =
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod,
		    FI->QuantityStorageDof_P->DefineQuantity->
		    IntegralQuantity.FMMIntegrationMethodIndex)) ->CriterionIndex ;
    
      
    List_Read(Problem_S.FMMGroup, FMMSrc, &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->A_L, i_Group, &Aggreg_M) ;
      List_Read(FMMmat_P->NumDof, i_Group, &NumDof_L) ;
      Nbr_DofList = List_Nbr(NumDof_L) ;

      if (Aggreg_M == NULL){		  
	Aggreg_M = (double**)Malloc(Nbr_DofList*sizeof(double*)) ;
	for (j = 0 ; j < Nbr_DofList ; j++)
	  Aggreg_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.Rsrc = 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 ;
	
	Current.Element = &Element ;
	Element.ElementSource = &Element ;
	Current.ElementSource = Element.ElementSource;
	
	QuantityStorageDof_P->NumLastElementForFunctionSpace = Element.Num ; 
	Get_DofOfElement
	  (&Element, QuantityStorageDof_P->FunctionSpace, QuantityStorageDof_P,
	   NULL) ;

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

	FI->xChangeOfCoordinatesDof =
	  (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->IntegralQuantityActive.Type_FormDof) ;
	  
	i = 0 ;

	while ((i <  List_Nbr(FI->IntegralQuantityActive.JacobianCase_L)) &&
	       ((j = ((struct JacobianCase *)List_Pointer(FI->IntegralQuantityActive.JacobianCase_L, i))
		 ->RegionIndex) >= 0) &&
	       !List_Search
	       (((struct Group *)List_Pointer(Problem_S.Group, j)) ->InitialList,
		&Element.ElementSource->Region, fcmp_int) )  i++ ;
    
	if (i == List_Nbr(FI->IntegralQuantityActive.JacobianCase_L))
	  Msg(GERROR, "Undefined Jacobian in Region %d", Element.Region);
	
	Element.JacobianCase = (struct JacobianCase*)List_Pointer(FI->IntegralQuantityActive.JacobianCase_L, 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, FMMIntegrationCaseDof_L, FMMCriterionIndexDof);
	  
	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,
			      FI->QuantityStorageDof_P->DefineQuantity->
			      IntegralQuantity.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.us, &Current.vs, &Current.ws, &weight) ;
	    
	    if (FI->Flag_ChangeCoord) {
	      Get_BFGeoElement(&Element, Current.us, Current.vs, Current.ws) ;
	      
	      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.xs = Current.ys = Current.zs = 0. ;
	    for (i = 0 ; i < Element.GeoElement->NbrNodes ; i++) {
	      Current.xs += Element.x[i] * Element.n[i] ;
	      Current.ys += Element.y[i] * Element.n[i] ;
	      Current.zs += Element.z[i] * Element.n[i] ;
	    }
	  

	    Factor = (FI->Flag_ChangeCoord) ? weight * fabs(Element.DetJac) : weight ;
	    
	    Cal_TermsforAggregation( QuantityStorageDof_P, Val0, &FunctionProd );
	    
	    Current.FMM.Flag_GF = FMM_AGGREGATION ;	
	    ((CASTF2V)QuantityStorageDof_P->DefineQuantity->IntegralQuantity.FunctionForFMM.Fct)
	      (&QuantityStorageDof_P->DefineQuantity->IntegralQuantity.FunctionForFMM,
	       GFValue, GFValue ) ;

	    for (i = 0 ; i < Nbr_Dof ; i++) {	    
	      if (QuantityStorageDof_P->BasisFunction[i].Dof->Type == DOF_UNKNOWN){
		xFunctionBFDof[i]
		  (&Element,
		   QuantityStorageDof_P->BasisFunction[i].NumEntityInElement+1,
		   Current.us, Current.vs, Current.ws, vBFuDof) ;

		vBFxDof[0] = vBFxDof[1]= vBFxDof[2] = 0. ;
		((void (*)(struct Element*, double*, double*))
		 FI->xChangeOfCoordinatesDof) (&Element, vBFuDof, vBFxDof) ;
				
		vBFxDofV.Type = Get_ValueFromForm(FI->Type_FormDof) ;
      		vBFxDofV.Val[0]         = vBFxDof[0] * Factor ;
		vBFxDofV.Val[1]         = vBFxDof[1] * Factor ;
		vBFxDofV.Val[2]         = vBFxDof[2] * Factor ;
		vBFxDofV.Val[MAX_DIM]   = 0. ;
		vBFxDofV.Val[MAX_DIM+1] = 0. ;
		vBFxDofV.Val[MAX_DIM+2] = 0. ;

		Apply_ConstantFactor(QuantityStorageDof_P, &vBFxDofV, &Val0) ;
		
		for (i_FMM = 0 ; i_FMM < NbrDir ; i_FMM++){ 		
		  ((CAST3V)FunctionProd)(&vBFxDofV, &GFValue[i_FMM], &BFGFValue[i_FMM]) ;
		}
		j = List_ISearch(NumDof_L, 
				 &QuantityStorageDof_P->BasisFunction[i].Dof->
				 Case.Unknown.NumDof, fcmp_int) ;

		if (j != -1) Cal_AddValueArray2DoubleArray(BFGFValue, Aggreg_M[j], NbrDir) ; 
		else  Msg(GERROR, "Wrong NumEqu %d for Aggreg",
			  QuantityStorageDof_P->BasisFunction[i].Dof->Case.Unknown.NumDof) ;
	      }/* if DOF_UNKNOWN */     
	    } /* for i Nbr_Dof */
	    
	  } /* for i_IntPoint ... */
	  break ; /* case GAUSS */
	    
	case ANALYTIC :
	  Current.FMM.Flag_GF = FMM_AGGREGATION ;
	  Current.FMM.Flag_Normal = 0;
	  if(((CASTF2V)QuantityStorageDof_P->DefineQuantity->IntegralQuantity.FunctionForFMM.Fct)
	     == ((CASTF2V)GF_NPxGradLaplace))	   
	    Current.FMM.Flag_Normal = 1 ;
	  
	  
	  Cal_ZeroValue(&Val0);
	  Cal_TermsforAggregation( QuantityStorageDof_P, Val0, &FunctionProd );
	  GF_FMMLaplacexForm( &Element, QuantityStorageDof_P, Nbr_Dof, 
			      (void (*)())xFunctionBFDof, NumDof_L, Val0, Aggreg_M );
	  break;
	}/*switch Integration_Case */	 
	
      } /* i_Element */
      List_Write(FMMmat_P->A_L, i_Group, &Aggreg_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 A g g r e g a t i o n                 */
/*-------------------------------------------------------------------------------*/

void Cal_TermsforAggregation( struct QuantityStorage * QuantityStorageDof_P,
			      struct Value Val0, void(**FunctionProd)()){

  struct WholeQuantity      * WholeQuantity_P0 ;


  GetDP_Begin("Cal_TermsForAggregation") ;


  WholeQuantity_P0 = (struct WholeQuantity*)
    List_Pointer(QuantityStorageDof_P->DefineQuantity->IntegralQuantity.WholeQuantity, 0) ;

  switch(QuantityStorageDof_P->DefineQuantity->IntegralQuantity.CanonicalWholeQuantity){
  case CWQ_GF :
    break ;
  case CWQ_GF_PSCA_DOF :
  case CWQ_GF_PVEC_DOF :
  case CWQ_DOF_PVEC_GF :	
    *FunctionProd = (WholeQuantity_P0+2)->Case.Operator.Function ;
    break ;
  case CWQ_GF_PSCA_EXP :
  case CWQ_GF_PVEC_EXP :
  case CWQ_EXP_PVEC_GF :
    *FunctionProd = (WholeQuantity_P0+2)->Case.Operator.Function ;
    Get_ValueOfExpression((struct Expression *)
			  List_Pointer(Problem_S.Expression,
				       QuantityStorageDof_P->DefineQuantity->IntegralQuantity.
				       ExpressionIndexForCanonical),
			  NULL, 0., 0., 0., &Val0) ;
    break ;
  case CWQ_EXP_TIME_GF_PSCA_DOF :
  case CWQ_EXP_TIME_GF_PVEC_DOF :
  case CWQ_EXP_PVEC_GF_PSCA_DOF :
  case CWQ_EXP_PVEC_GF_PVEC_DOF :
    *FunctionProd = (WholeQuantity_P0+4)->Case.Operator.Function ;
    Get_ValueOfExpression((struct Expression *)
			  List_Pointer(Problem_S.Expression,
				       QuantityStorageDof_P->DefineQuantity->IntegralQuantity.
				       ExpressionIndexForCanonical),
			  NULL, 0., 0., 0., &Val0) ;
    break;
  case CWQ_FCT_TIME_GF_PSCA_DOF :
  case CWQ_FCT_TIME_GF_PVEC_DOF :
  case CWQ_FCT_PVEC_GF_PSCA_DOF :
  case CWQ_FCT_PVEC_GF_PVEC_DOF :
    *FunctionProd = (WholeQuantity_P0+4)->Case.Operator.Function ;
    ((CASTF2V)QuantityStorageDof_P->DefineQuantity->IntegralQuantity.AnyFunction.Fct)
      (&QuantityStorageDof_P->DefineQuantity->IntegralQuantity.AnyFunction, &Val0, &Val0) ;	 
    break;
  default: Msg(GERROR, "Unknown type of Integral Quantity for FMM");
  }
  
  
  GetDP_End ;
  
} 



#undef CAST3V 
#undef CASTF2V