fmm_groups.c

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

      FoundHome = 0;
      
      for (j = 0 ; j < NbrDivX ; j++){
	for (k = 0 ; k < NbrDivY ; k++){
	  for (l = 0 ; l < NbrDivZ ; l++){
	    x1 = -tol+Xmin +(double)j*(Xmax-Xmin)/(double)NbrDivX;
	    x2 =  tol+Xmin +(double)(j+1)*(Xmax-Xmin)/(double)NbrDivX;
	    y1 = -tol+Ymin +(double)k*(Ymax-Ymin)/(double)NbrDivY;
	    y2 =  tol+Ymin +(double)(k+1)*(Ymax-Ymin)/(double)NbrDivY;
	    z1 = -tol+Zmin +(double)l*(Zmax-Zmin)/(double)NbrDivZ;
	    z2 =  tol+Zmin +(double)(l+1)*(Zmax-Zmin)/(double)NbrDivZ;

	    if( Current.Region == (Elements_P0+iElm)->Region ){
	      Geo_SetCentroidCoordinates((Elements_P0+iElm)->Num, Centroid);
	      
	      if (x1 <= Centroid[0] && Centroid[0] <= x2 &&
		  y1 <= Centroid[1] && Centroid[1] <= y2 &&
		  z1 <= Centroid[2] && Centroid[2] <= z2 ){
		if(!NbrElmsDiv[j][k][l]){
		  NbrElmsDiv[j][k][l] = 1;
		  ElmsDiv[j][k][l][0] = iElm ;
		} else {
		  NbrElmsDiv[j][k][l] += 1;
		  ElmsDiv[j][k][l] = (int *)realloc(ElmsDiv[j][k][l], NbrElmsDiv[j][k][l]*sizeof(int));
		  ElmsDiv[j][k][l][NbrElmsDiv[j][k][l]-1] = iElm;	    
		}
		FoundHome = 1; break;
	      }
	    }/* if Region */
	    if(FoundHome)break;
	  }
	  if(FoundHome)break;
	}
	if(FoundHome)break;
      }
      if(!FoundHome && (Current.Region == (Elements_P0+iElm)->Region)) 
	Msg(WARNING,"Element %d has not been assigned to a FMMGroup", iElm);
    }/* Elms NbrElms */ 

    NbrFMMGroupsInRegion = 0;
    for (j = 0 ; j < NbrDivX ; j++)
      for (k = 0 ; k < NbrDivY ; k++)
	for (l = 0 ; l < NbrDivZ ; l++)
	  if(NbrElmsDiv[j][k][l]){
	    Geo_InitFMMData(&FMMData_S) ;
	    FMMData_S.Element = List_Create(NbrElmsDiv[j][k][l], 1, sizeof(int)) ;
	    for(m = 0 ; m < NbrElmsDiv[j][k][l] ; m++){
	      List_Add(FMMData_S.Element, &(Elements_P0 + ElmsDiv[j][k][l][m])->Num) ;
	      Geo_SetCentroidCoordinates((Elements_P0+ElmsDiv[j][k][l][m])->Num, Centroid) ;

	      FMMData_S.Xgc += Centroid[0] ;
	      FMMData_S.Ygc += Centroid[1] ;
	      FMMData_S.Zgc += Centroid[2] ;
	      
	      (Elements_P0 + ElmsDiv[j][k][l][m])->FMMGroup = iFMMGroup ;
	    }	
	    
	    FMMData_S.Group = iFMMGroup;  
	    FMMData_S.Xgc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Ygc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Zgc /= NbrElmsDiv[j][k][l];
	    List_Add(FMMGroup_L, &FMMData_S) ;
	    
	    iFMMGroup +=1 ;
	    NbrFMMGroupsInRegion += 1 ;
	  }/* if( NbrElmsDiv[j][k][l] ) */    
  
    
    Free(NbrElmsDiv);
    Free(ElmsDiv);

    Msg(INFO,"SupportIndex = %d  Region = %d  NbrElmsInRegion = %d  NbrFMMGroupsInRegion = %d",
	InSupport, Current.Region, NbrElmsInRegion, NbrFMMGroupsInRegion) ; 
  }/* iGroup NbrGroupInSuppport*/
 
  Msg(INFO,"SupportIndex = %d  NbrElmsInSupport = %d  NbrFMMGroupsInSupport = %d",
      InSupport, NbrElmsInSupport, iFMMGroup) ; 

  FMMGroup_S.InIndex = InSupport ;
  FMMGroup_S.List = FMMGroup_L ;
  List_Add(Problem_S.FMMGroup, &FMMGroup_S ); 


  GetDP_End ;
}


#if 0 /* FIXME: this is just a sketch -- is not meant to work or even compile */
  
void  Geo_CreateMultilevelFMMGroup( int InSupport, struct GeoData *GeoData_P, double k0 ){

  int NbrGroupsInSupport, iGroup, NbrFMMGroupsInRegion ;
  int NbrElms, NbrElmsInRegion, NbrElmsInSupport ;
  int j, k, l, m, FoundHome, iFMMGroup, iElm, iFMMParentGroup ;
  int ***NbrElmsDiv, ****ElmsDiv, ***NbrGrsDiv, ****GrsDiv ;
  double  Xmax, Xmin, Ymax, Ymin, Zmax, Zmin ;
  int NbrDivX, NbrDivY, NbrDivZ;
  int NbrDivX0 = 1, NbrDivY0 = 1, NbrDivZ0 = 1, NbrDivX1 = 1, NbrDivY1 = 1, NbrDivZ1 = 1  ;
  double  x1, x2, y1, y2, z1, z2, Centroid[3], tol ;
  struct Value DivXYZ ;
  struct Geo_Element *Elements_P0 ;
  struct FMMData FMMData_S ;
  struct FMMData FMMParentData_S ;
  struct FMMGroup FMMGroup_S ;

  List_T * FMMGroup_L ;  
  List_T * FMMParentGroup_L ;  


  GetDP_Begin("Geo_CreateMultilevelFMMGroup");

  tol = 1e-6;

  NbrElms = List_Nbr(GeoData_P->Elements) ;
  Elements_P0 =  (struct Geo_Element*)List_Pointer(GeoData_P->Elements, 0) ;

  FMMGroup_L = List_Create (10, 10, sizeof(struct FMMData)) ; 
  FMMParentGroup_L = List_Create (10, 10, sizeof(struct FMMData)) ;   

  NbrElmsInSupport = 0 ;  
  NbrGroupsInSupport = List_Nbr(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList ); 

  iFMMGroup = 0 ;
  iFMMParentGroup = 0 ;
  
  for (iGroup = 0 ; iGroup < NbrGroupsInSupport ; iGroup++){
    List_Read(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList, iGroup, &Current.Region) ;    
    Get_ValueOfExpressionByIndex(Current.FMM.DivXYZIndex, NULL, 0., 0., 0., &DivXYZ) ;
    NbrDivX0 = (int)DivXYZ.Val[0] ;
    NbrDivY0 = (int)DivXYZ.Val[1] ;
    NbrDivZ0 = (int)DivXYZ.Val[2] ;

    NbrDivX1 = (int)DivXYZ.Val[MAX_DIM] ;
    NbrDivY1 = (int)DivXYZ.Val[MAX_DIM+1] ;
    NbrDivZ1 = (int)DivXYZ.Val[MAX_DIM+2] ;

    Msg(INFO,"Div. Parents -> Region = %d DivX0 = %d  DivY0 = %d  DivZ0 = %d", Current.Region, NbrDivX0, NbrDivY0, NbrDivZ0) ; 
    Msg(INFO,"Div. Children -> DivX1 = %d  DivY1 = %d  DivZ1 = %d", NbrDivX1, NbrDivY1, NbrDivZ1) ; 
    
    NbrElmsInRegion = 0 ; 
    Xmax = Ymax = Zmax = Xmin = Ymin = Zmin = 0. ;

    for(iElm = 0 ; iElm < NbrElms ; iElm++){
      if((Elements_P0+iElm)->Region == Current.Region){
	NbrElmsInRegion++ ; 
	Geo_SetXYZMinMaxInElement((Elements_P0+iElm)->Num, &Xmin, &Ymin, &Zmin, &Xmax, &Ymax, &Zmax ) ;      
      }      
    }

    NbrElmsInSupport += NbrElmsInRegion ; 

    NbrElmsDiv = (int ***)Malloc(NbrDivX*sizeof(int**));
    for (j = 0 ; j < NbrDivX ; j++){
      NbrElmsDiv[j] = (int **)Calloc(NbrDivY, sizeof(int*));
      for (k = 0 ; k < NbrDivY ; k++)
	NbrElmsDiv[j][k] = (int *)Calloc(NbrDivZ, sizeof(int));
    }
    
    ElmsDiv = (int ****)Malloc(NbrDivX*sizeof(int***));
    for (j = 0 ; j < NbrDivX ; j++){
      ElmsDiv[j] = (int ***)Malloc(NbrDivY*sizeof(int**)); 
      for (k = 0 ; k < NbrDivY ; k++){
	ElmsDiv[j][k] = (int **)Malloc(NbrDivZ*sizeof(int*));
	for (l = 0 ; l < NbrDivZ ; l++)
	  ElmsDiv[j][k][l] = (int *)Malloc(sizeof(int));
      }
    }


    NbrGrsDiv = (int ***)Malloc(NbrDivX0*sizeof(int**));
    for (j = 0 ; j < NbrDivX0 ; j++){
      NbrGrsDiv[j] = (int **)Calloc(NbrDivY0, sizeof(int*));
      for (k = 0 ; k < NbrDivY0 ; k++)
	NbrGrsDiv[j][k] = (int *)Calloc(NbrDivZ0, sizeof(int));
    }


    GrsDiv = (int ****)Malloc(NbrDivX0*sizeof(int***));
    for (j = 0 ; j < NbrDivX0 ; j++){
      GrsDiv[j] = (int ***)Malloc(NbrDivY0*sizeof(int**)); 
      for (k = 0 ; k < NbrDivY0 ; k++){
	GrsDiv[j][k] = (int **)Malloc(NbrDivZ0*sizeof(int*));
	for (l = 0 ; l < NbrDivZ0 ; l++)
	  GrsDiv[j][k][l] = (int *)Malloc(sizeof(int));
      }
    }


    NbrDivX = NbrDivX0 * NbrDivX1 ;
    NbrDivY = NbrDivY0 * NbrDivY1 ;
    NbrDivZ = NbrDivZ0 * NbrDivZ1 ;

    for (iElm = 0 ; iElm < NbrElms ; iElm++){ 
      FoundHome = 0;
      
      for (j = 0 ; j < NbrDivX ; j++){
	for (k = 0 ; k < NbrDivY ; k++){
	  for (l = 0 ; l < NbrDivZ ; l++){
	    x1 = -tol+Xmin +(double)j*(Xmax-Xmin)/(double)NbrDivX;
	    x2 =  tol+Xmin +(double)(j+1)*(Xmax-Xmin)/(double)NbrDivX;
	    y1 = -tol+Ymin +(double)k*(Ymax-Ymin)/(double)NbrDivY;
	    y2 =  tol+Ymin +(double)(k+1)*(Ymax-Ymin)/(double)NbrDivY;
	    z1 = -tol+Zmin +(double)l*(Zmax-Zmin)/(double)NbrDivZ;
	    z2 =  tol+Zmin +(double)(l+1)*(Zmax-Zmin)/(double)NbrDivZ;

	    if( Current.Region == (Elements_P0+iElm)->Region ){
	      Geo_SetCentroidCoordinates((Elements_P0+iElm)->Num, Centroid);
	      
	      if (x1 <= Centroid[0] && Centroid[0] <= x2 &&
		  y1 <= Centroid[1] && Centroid[1] <= y2 &&
		  z1 <= Centroid[2] && Centroid[2] <= z2 ){
		if(!NbrElmsDiv[j][k][l]){
		  NbrElmsDiv[j][k][l] = 1;
		  ElmsDiv[j][k][l][0] = iElm ;
		} else {
		  NbrElmsDiv[j][k][l] += 1;
		  ElmsDiv[j][k][l] = (int *)realloc(ElmsDiv[j][k][l], NbrElmsDiv[j][k][l]*sizeof(int));
		  ElmsDiv[j][k][l][NbrElmsDiv[j][k][l]-1] = iElm;	    
		}
		FoundHome = 1; break;
	      }
	    }/* if Region */
	    if(FoundHome)break;
	  }
	  if(FoundHome)break;
	}
	if(FoundHome)break;
      }
      if(!FoundHome && (Current.Region == (Elements_P0+iElm)->Region)) 
	Msg(WARNING,"Element %d has not been assigned to a FMMGroup", iElm);
    }/* Elms NbrElms */ 


    NbrFMMGroupsInRegion = 0;
    for (j = 0 ; j < NbrDivX ; j++)
      for (k = 0 ; k < NbrDivY ; k++)
	for (l = 0 ; l < NbrDivZ ; l++)
	  if(NbrElmsDiv[j][k][l]){
	    Geo_InitFMMData(&FMMData_S) ;
	    FMMData_S.Element = List_Create(NbrElmsDiv[j][k][l], 1, sizeof(int)) ;
	    for(m = 0 ; m < NbrElmsDiv[j][k][l] ; m++){
	      List_Add(FMMData_S.Element, &(Elements_P0 + ElmsDiv[j][k][l][m])->Num) ;
	      Geo_SetCentroidCoordinates((Elements_P0+ElmsDiv[j][k][l][m])->Num, Centroid) ;

	      FMMData_S.Xgc += Centroid[0] ;
	      FMMData_S.Ygc += Centroid[1] ;
	      FMMData_S.Zgc += Centroid[2] ;
	      
	      (Elements_P0 + ElmsDiv[j][k][l][m])->FMMGroup = iFMMGroup ;
	    }	
	    
	    FMMData_S.Group = iFMMGroup;  
	    FMMData_S.Xgc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Ygc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Zgc /= NbrElmsDiv[j][k][l];

	    List_Add(FMMGroup_L, &FMMData_S) ;
	    
	    iFMMGroup +=1 ;
	    NbrFMMGroupsInRegion += 1 ;
	  }/* if( NbrElmsDiv[j][k][l] ) */ 
  
    Free(NbrElmsDiv);
    Free(ElmsDiv); 

    Msg(INFO,"SupportIndex = %d  Region = %d  NbrElmsInRegion = %d  NbrFMMGroupsInRegion = %d",
	InSupport, Current.Region, NbrElmsInRegion, NbrFMMGroupsInRegion) ;
  }/* iGroup NbrGroupInSuppport*/
  
  Msg(INFO,"SupportIndex = %d  NbrElmsInSupport = %d  NbrFMMGroupsInSupport = %d",
      InSupport, NbrElmsInSupport, iFMMGroup) ; 

  for (iFMMGroup = 0 ; iFMMGroup < List_Nbr(FMMGroup_L) ; iFMMGroup++){       
    FoundHome = 0;
    List_Read(FMMGroup_L, iFMMGroup, &FMMData_S); 
    for (j = 0 ; j < NbrDivX0 ; j++){
      for (k = 0 ; k < NbrDivY0 ; k++){
	for (l = 0 ; l < NbrDivZ0 ; l++){
	  x1 = -tol+Xmin +(double)j*(Xmax-Xmin)/(double)NbrDivX0;
	  x2 =  tol+Xmin +(double)(j+1)*(Xmax-Xmin)/(double)NbrDivX0;
	  y1 = -tol+Ymin +(double)k*(Ymax-Ymin)/(double)NbrDivY0;
	  y2 =  tol+Ymin +(double)(k+1)*(Ymax-Ymin)/(double)NbrDivY0;
	  z1 = -tol+Zmin +(double)l*(Zmax-Zmin)/(double)NbrDivZ0;
	  z2 =  tol+Zmin +(double)(l+1)*(Zmax-Zmin)/(double)NbrDivZ0;
	  
	  if (x1 <=  FMMData_S.Xgc && FMMData_S.Xgc <= x2 &&
	      y1 <=  FMMData_S.Ygc && FMMData_S.Ygc <= y2 &&
	      z1 <=  FMMData_S.Zgc && FMMData_S.Zgc <= z2 ){
	    if(!NbrGrsDiv[j][k][l]){
	      NbrGrsDiv[j][k][l] = 1;
	      GrsDiv[j][k][l][0] = iFMMGroup ;
	    }
	    else {
	      NbrGrsDiv[j][k][l] += 1;
	      GrsDiv[j][k][l] = (int *)realloc(GrsDiv[j][k][l], NbrGrsDiv[j][k][l]*sizeof(int));
	      GrsDiv[j][k][l][NbrGrsDiv[j][k][l]-1] = iFMMGroup;
	    }
	    FoundHome = 1; break;
	  }	    
	  if(FoundHome)break;
	}
	if(FoundHome)break;
      }
      if(FoundHome)break;
    }
  } /* Grs in List_Nbr(FMMGroup_L) */ 
  
  for (j = 0 ; j < NbrDivX0 ; j++)
    for (k = 0 ; k < NbrDivY0 ; k++)
      for (l = 0 ; l < NbrDivZ0 ; l++)
	if(NbrGrsDiv[j][k][l]){
	  Geo_InitFMMData(&FMMParentData_S) ;
	  FMMParentData_S.Element = List_Create(NbrGrsDiv[j][k][l], 1, sizeof(struct FMMData)) ;
	  for(m = 0 ; m < NbrGrsDiv[j][k][l] ; m++){
	    List_Read(FMMGroup_L, GrsDiv[j][k][l][m], &FMMData_S); 
	    List_Add(FMMParentData_S.Element, &FMMData_S) ;
	    
	    FMMParentData_S.Xgc +=  FMMData_S.Xgc ;
	    FMMParentData_S.Ygc +=  FMMData_S.Ygc ;
	    FMMParentData_S.Zgc +=  FMMData_S.Zgc ;   
	  }	
	  
	  FMMParentData_S.Group = iFMMParentGroup;  
	  FMMParentData_S.Xgc /= NbrGrsDiv[j][k][l];
	  FMMParentData_S.Ygc /= NbrGrsDiv[j][k][l];
	  FMMParentData_S.Zgc /= NbrGrsDiv[j][k][l];
	  
	  List_Add(FMMParentGroup_L, &FMMParentData_S) ;
	  
	  iFMMParentGroup +=1 ;
	}/* if( NbrGrsDiv[j][k][l] ) */
  
  
  Free(NbrGrsDiv);
  Free(GrsDiv);
  
  List_Reset(FMMGroup_L);
  
  FMMGroup_S.InIndex = InSupport ;
  FMMGroup_S.List = FMMParentGroup_L ;
  List_Add(Problem_S.FMMGroup, &FMMGroup_S ); 
  
  GetDP_End ;
}

#endif


void Get_InFMMGroupList( int Index_Formulation, struct GeoData *GeoData_P ){

  int  i_EquTerm, Nbr_EquationTerm, i_Observation, i_Source, FlagError, DIM ;
  double k0 = -1. ;
  struct Formulation    * Formulation_P ;
  struct EquationTerm   * EquationTerm_P0, * EquationTerm_P ;
  struct DefineQuantity * DefineQuantityDof_P, * DefineQuantity_P0 ;
  struct FMMmat         * FMMmat_P, FMMmat_S ;
  List_T *InIndex ;

  GetDP_Begin("Get_InFMMGroupList") ;

  Formulation_P = (struct Formulation*)List_Pointer(Problem_S.Formulation, Index_Formulation) ;
  EquationTerm_P0 = (struct EquationTerm*)List_Pointer(Formulation_P->Equation, 0) ;
  Nbr_EquationTerm = List_Nbr(Formulation_P->Equation) ;
  DefineQuantity_P0 = (struct DefineQuantity*)List_Pointer(Formulation_P->DefineQuantity, 0) ;

  InIndex = List_Create( 1, 1, sizeof(int)) ;

  Problem_S.FMMGroup = List_Create (1, 1, sizeof(struct FMMGroup)) ;

  Current.DofData->FMM_Matrix = NULL ;

  /* Current.FMM.Dimension = DIM = (int)((GeoData_P)->Dimension) ; */
  
  for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
    EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
    EquationTerm_P->Case.LocalTerm.FMMObservation =
      EquationTerm_P->Case.LocalTerm.FMMSource = EquationTerm_P->Case.LocalTerm.iFMMEqu = -1 ;

    if ( EquationTerm_P->Type == GALERKIN && EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex > 1 ) {
      DefineQuantityDof_P = DefineQuantity_P0 + EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof ;
      if( DefineQuantityDof_P->Type==INTEGRALQUANTITY){

	if(DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.NbrParameters == 2)
	  k0 = DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Para[1] ;/* Helmoltz Case */
	  DIM = (int)DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Para[0] ;
	  Current.FMM.Dimension = DIM ;

	if( List_PQuery(InIndex, &EquationTerm_P->Case.LocalTerm.InIndex, fcmp_int) == NULL ){
	  List_Add( InIndex, &EquationTerm_P->Case.LocalTerm.InIndex ) ;
	  Geo_CreateFMMGroup( EquationTerm_P->Case.LocalTerm.InIndex, GeoData_P, k0) ;
	}

	if( List_PQuery(InIndex, &DefineQuantityDof_P->IntegralQuantity.InIndex, fcmp_int) == NULL ){
	  List_Add( InIndex, &DefineQuantityDof_P->IntegralQuantity.InIndex ) ;
	  Geo_CreateFMMGroup( DefineQuantityDof_P->IntegralQuantity.InIndex, GeoData_P, k0 ) ;
	}

	

	if (Current.DofData->FMM_Matrix == NULL)
	  Current.DofData->FMM_Matrix = List_Create(1, 1, sizeof(struct FMMmat)) ;
	
	i_Source = EquationTerm_P->Case.LocalTerm.FMMSource = 
	  List_ISearch(InIndex, &DefineQuantityDof_P->IntegralQuantity.InIndex, fcmp_int) ;
	i_Observation = EquationTerm_P->Case.LocalTerm.FMMObservation = 
	  List_ISearch(InIndex, &EquationTerm_P->Case.LocalTerm.InIndex, fcmp_int) ;

	InitFMMmatrix(i_Source, i_Observation, &FMMmat_S) ;
	List_Add( Current.DofData->FMM_Matrix, &FMMmat_S) ;
 	EquationTerm_P->Case.LocalTerm.iFMMEqu = List_Nbr(Current.DofData->FMM_Matrix)-1 ;
	
	Get_GroupNeighbours( EquationTerm_P->Case.LocalTerm.iFMMEqu) ;

	FMMmat_P = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, EquationTerm_P->Case.LocalTerm.iFMMEqu) ;

	if (DIM == _2D)
	  Get_FunctionForFunction( FMMProd_Function2D, DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Fct,
				   &FlagError, &(FMMmat_P->FctProd) ) ;
	else
	  Get_FunctionForFunction( FMMProd_Function3D, DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Fct,
				   &FlagError, &(FMMmat_P->FctProd) ) ;

	FMMmat_P->GFx = &(DefineQuantityDof_P->IntegralQuantity.FunctionForFMM) ;

	
      }/* Quantity */
    }/* If Galerkin */
  }/* Equation_Term */

    
  Init_CurrentFMMData ( k0 ) ;
  List_Delete(InIndex);

  GetDP_End ; 

}