fmm_groups.c

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

      }
      List_Reset(FMMmat_P->A_L);
      List_Reset(FMMmat_P->D_L);
    }

  }
  
  for (i_FMMEqu = 0 ; i_FMMEqu < NbrFMMEqu ; i_FMMEqu++){
    FMMmat_P  = FMMmat_P0 + i_FMMEqu ;  
    i_Obs = FMMmat_P->Obs ;
    i_Src = FMMmat_P->Src ;
    i = 0 ; while (i < i_FMMEqu) { if ((FMMmat_P0+i)->Obs == i_Obs  &&  (FMMmat_P0+i)->Src == i_Src) break ;i++;}

    if (i < i_FMMEqu ){ 
      FMMmat_P->NG_L = (FMMmat_P0+i)->NG_L ;   
      FMMmat_P->FG_L = (FMMmat_P0+i)->FG_L ;
      FMMmat_P->Nd_L = (FMMmat_P0+i)->Nd_L ;

      List_Read(Problem_S.FMMGroup, i_Src, &FMMGroup_S) ;
      NbrGroupSrc =  NbrGroupObs = List_Nbr(FMMGroup_S.List) ;
      if (i_Src != i_Obs){
	List_Read(Problem_S.FMMGroup, i_Obs, &FMMGroup_S) ;
	NbrGroupObs = List_Nbr(FMMGroup_S.List) ;
      }   
      
      if(Flag_FMMDA){
	for(i = 0 ; i < NbrGroupSrc ; i++) List_Add(FMMmat_P->A_L, &Mat) ;
	for(i = 0 ; i < NbrGroupObs ; i++) List_Add(FMMmat_P->D_L, &Mat) ;
      }
      
    }
    else{
      List_Read(Problem_S.FMMGroup, i_Src, &FMMGroup_S) ;
      NbrGroupSrc =  NbrGroupObs = List_Nbr(FMMGroup_S.List) ;
      FMMDataGSrc_P0 = FMMDataGObs_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;
      
      if (i_Src != i_Obs){
	List_Read(Problem_S.FMMGroup, i_Obs, &FMMGroup_S) ;
	NbrGroupObs = List_Nbr(FMMGroup_S.List) ;
	FMMDataGObs_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;
      }

      
      if(Flag_FMMDA){
      for(i = 0 ; i < NbrGroupSrc ; i++) List_Add(FMMmat_P->A_L, &Mat) ;
      for(i = 0 ; i < NbrGroupObs ; i++) List_Add(FMMmat_P->D_L, &Mat) ;
      }
      

      for ( i = 0 ; i < NbrGroupSrc; i++){
	XSrc = (FMMDataGSrc_P0+i)->Xgc ;
	YSrc = (FMMDataGSrc_P0+i)->Ygc ;
	ZSrc = (FMMDataGSrc_P0+i)->Zgc ;
	
	NGi = List_Create(5,  2, sizeof(int)) ;
	FGi = List_Create(10, 2, sizeof(int)) ;
	Ndi = List_Create(10, 2, sizeof(int)) ;
	/* Dfar = Current.FMM.far * (FMMDataGSrc_P0+i)->Rmax *2 ; */
	
	for ( j = 0 ; j < NbrGroupObs; j++){
	  d = sqrt(SQU((FMMDataGObs_P0+j)->Xgc-XSrc) + SQU((FMMDataGObs_P0+j)->Ygc-YSrc) +
		   SQU((FMMDataGObs_P0+j)->Zgc-ZSrc)) ;
	  
	  if(d < Dfar) List_Add(NGi, &j) ;
	  else {
	    List_Add(FGi, &j) ;

	    /* Ndir = FMM_SetTruncation((FMMDataGSrc_P0+i)->Rmax,(FMMDataGObs_P0+j)->Rmax,
	                                d, Current.FMM.Dimension) ; */
	    Ndir = FMM_SetTruncationOLD((FMMDataGSrc_P0+i)->Rmax, d, Current.FMM.Dimension) ;
	    List_Add(Ndi, &Ndir) ;
	  }
	}/* NbrGroupObs */

	List_Add( FMMmat_P->NG_L, &NGi );
	List_Add( FMMmat_P->FG_L, &FGi );
	List_Add( FMMmat_P->Nd_L, &Ndi );
      }/* NbrGroupSrc */
    }
  }/* i_FMMEqu */

  GetDP_End;

}


void  Geo_SetXYZMinMaxInElement(int Num, double *Xmin, double *Ymin, double *Zmin,
				double *Xmax, double *Ymax, double *Zmax){
 
  int j ;
  double x[NBR_MAX_NODES_IN_ELEMENT], y[NBR_MAX_NODES_IN_ELEMENT], z[NBR_MAX_NODES_IN_ELEMENT];

  struct Geo_Element *Element ;  

  GetDP_Begin("Geo_SetXYZMinMaxInElement");

  Element =  Geo_GetGeoElementOfNum(Num) ;    
  Geo_GetNodesCoordinates(Element->NbrNodes, Element->NumNodes,
			  x, y, z);

  for (j = 0;  j < Element->NbrNodes ; j++){
    *Xmax = (*Xmax >= x[j]) ? *Xmax : x[j] ;
    *Ymax = (*Ymax >= y[j]) ? *Ymax : y[j] ;
    *Zmax = (*Zmax >= z[j]) ? *Zmax : z[j] ;
    *Xmin = (*Xmin <= x[j]) ? *Xmin : x[j] ;
    *Ymin = (*Ymin <= y[j]) ? *Ymin : y[j] ;
    *Zmin = (*Zmin <= z[j]) ? *Zmin : z[j] ;
  }
  
  GetDP_End ;
}



double Geo_SetMaxEdgeLength(int Num){
 
  int j ;
  double D = 0. ,Dmax =0. ;
  double x[NBR_MAX_NODES_IN_ELEMENT], y[NBR_MAX_NODES_IN_ELEMENT], z[NBR_MAX_NODES_IN_ELEMENT];

  struct Geo_Element *Element ;  


  GetDP_Begin("Geo_SetMaxEdgeLength");

  Element =  Geo_GetGeoElementOfNum(Num) ;    
  Geo_GetNodesCoordinates(Element->NbrNodes, Element->NumNodes,
			  x, y, z);
  
  for (j = 0;  j < Element->NbrNodes ; j++){
    if (j+1<Element->NbrNodes)
      D = sqrt(SQU(x[j]-x[j+1])+ SQU(y[j]-y[j+1])+ SQU( z[j]-z[j+1])) ;
  
    Dmax = ( D > Dmax ) ? D : Dmax ;
  }
  
  GetDP_Return(Dmax) ;
}


double Geo_GetMaxDistance2Elms(int Numi, int Numj){
 
  int i, j ;
  double D = 0. ,Dmax =0. ;
  double x1[NBR_MAX_NODES_IN_ELEMENT], y1[NBR_MAX_NODES_IN_ELEMENT], z1[NBR_MAX_NODES_IN_ELEMENT];
  double x2[NBR_MAX_NODES_IN_ELEMENT], y2[NBR_MAX_NODES_IN_ELEMENT], z2[NBR_MAX_NODES_IN_ELEMENT];

  struct Geo_Element *Elementi, *Elementj ;  

  GetDP_Begin("Geo_GetMaxDistance2Elms");

  Elementi =  Geo_GetGeoElementOfNum(Numi) ;    
  Geo_GetNodesCoordinates(Elementi->NbrNodes, Elementi->NumNodes,
			  x1, y1, z1);

  Elementj =  Geo_GetGeoElementOfNum(Numj) ;    
  Geo_GetNodesCoordinates(Elementj->NbrNodes, Elementj->NumNodes,
			  x2, y2, z2);
  
  for (i = 0;  i < Elementi->NbrNodes ; i++)
    for (j = 0;  j < Elementj->NbrNodes ; j++){
      D = sqrt(SQU(x1[i]-x2[j])+ SQU(y1[i]-y2[j])+ SQU(z1[i]-z2[j])) ;  
      Dmax = ( D > Dmax ) ? D : Dmax ;
    }
  
  GetDP_Return(Dmax) ;
}



int FMM_NbrSpectralDirections( double k0 ){

  int i, j, i_Support, NbrInSupport, NbrGroups, NbrElmsGroupi ;
  int NumElmj, NumElmjp ;
  double D = 0., Daux = 0., Centroidj[3], Centroidjp[3] ; 
  struct FMMData  *FMMData_P0 ;
  struct FMMGroup  FMMGroup_S ;

  GetDP_Begin("FMM_NbrSpectralDirections");

  NbrInSupport = List_Nbr(Problem_S.FMMGroup) ;

  for( i_Support = 0 ; i_Support < NbrInSupport ; i_Support++ ){
    List_Read(Problem_S.FMMGroup, i_Support, &FMMGroup_S) ;
    NbrGroups = List_Nbr(FMMGroup_S.List) ;
    FMMData_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;

    for(i = 0; i < NbrGroups; i++ ){
      NbrElmsGroupi = List_Nbr((FMMData_P0+i)->Element) ;
      
      if (NbrElmsGroupi==1){ 
	List_Read((FMMData_P0+i)->Element, 0, &NumElmj);
	Daux = Geo_SetMaxEdgeLength(NumElmj);
      }
      else{
	for (j = 0 ; j < NbrElmsGroupi ; j++){
	  if ((j+1) < NbrElmsGroupi){
	    List_Read((FMMData_P0+i)->Element, j, &NumElmj);
	    List_Read((FMMData_P0+i)->Element, j+1, &NumElmjp);
	    
	    Geo_SetCentroidCoordinates(NumElmj, Centroidj);
	    Geo_SetCentroidCoordinates(NumElmjp, Centroidjp);
	    
	    Daux = sqrt(SQU(Centroidj[0]-Centroidjp[0])+ SQU(Centroidj[1]-Centroidjp[1])+ SQU(Centroidj[2]-Centroidjp[2])) ;       
	  }
	  D = (Daux>D) ? Daux : D ;
	}
      }
      D = (Daux>D) ? Daux : D ;
    }    
  }
  
  Msg(INFO,"D_FMM = %.8f", D);
  
  GetDP_Return( (int)ceil( k0*D+5.*log(k0*D+PI) ) );
}




int FMM_SetTruncation( double Rsrc, double Robs, double D, int Dimension){

  int NbrDir=0, i, j, k, exp ;
  double Precision ; 

  GetDP_Begin("FMM_SetTruncation");

  Precision = Current.FMM.Precision ;
  exp = - (int)floor(log10(Precision)) ;

  i = (int)floor(Rsrc/D/0.05)-1 ; 
  j = (int)floor(Robs/D/0.05)-1 ;
  if (i<0) i = 0 ;
  if (j<0) j = 0 ;
  k = ((i + j * 7)>48)?48: i+j*7 ;


  if (Dimension == _2D){
    switch(exp){
    case 3 :  NbrDir = t2D1e_3[k] ; break ;  
    case 4 :  NbrDir = t2D1e_4[k] ; break ;  
    case 5 :  NbrDir = t2D1e_5[k] ; break ;  
    case 6 :  NbrDir = t2D1e_6[k] ; break ;  
    case 7 :  NbrDir = t2D1e_7[k] ; break ;  
    case 8 :  NbrDir = t2D1e_8[k] ; break ;  
    case 9 :  NbrDir = t2D1e_9[k] ; break ;
    default: 
      Msg(GERROR, "Precision for FMM truncation can only be: 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9");
    }
    if (NbrDir > 11){ if (Flag_FMM == 1) Msg(WARNING,"NbrDir = %d changed to 11 for memory requirements", NbrDir) ; NbrDir = 11 ;}
  }
  else{
    switch(exp){
    case 3 :  NbrDir = t3D1e_3[k] ; break ;  
    case 4 :  NbrDir = t3D1e_4[k] ; break ;  
    case 5 :  NbrDir = t3D1e_5[k] ; break ;  
    case 6 :  NbrDir = t3D1e_6[k] ; break ;  
    case 7 :  NbrDir = t3D1e_7[k] ; break ;  
    case 8 :  NbrDir = t3D1e_8[k] ; break ;  
    case 9 :  NbrDir = t3D1e_9[k] ; break ;
    default: 
      Msg(GERROR, "Precision for FMM truncation can only be: 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9");
  }
    if (NbrDir > 15){ 
      if (Flag_FMM == 1) 
	Msg(WARNING,"NbrDir = %d changed to 15 for memory requirements Rsrc/D %g Robs/D %g", NbrDir,Rsrc/D,Robs/D) ; NbrDir = 15 ;}
  }
 

  if (Flag_FMM == 2)  
    NbrDir = (NbrDir > Current.FMM.NbrDir) ? Current.FMM.NbrDir : NbrDir ; 
  else
    Current.FMM.NbrDir = (NbrDir>Current.FMM.NbrDir) ? NbrDir : Current.FMM.NbrDir ;
  
  GetDP_Return(NbrDir);
}


int FMM_SetTruncationOLD( double Rmax, double D, int Dimension){

  int NbrDir ;
  double Precision ; 

  GetDP_Begin("FMM_SetTruncation");

 /* Rmax == Radius of the circle/sphere enclosing the FMM group */  
  Precision = Current.FMM.Precision ;

  if (Dimension == _2D){
    NbrDir =  (int)floor(-log(Precision)/log(D/Rmax)) ; /* Truncation parameter */
    if (NbrDir > 15){ if (Flag_FMM == 1) Msg(WARNING,"NbrDir = %d changed to 15 for memory requirements", NbrDir) ; NbrDir = 15 ;}
  }
  else{
    NbrDir =  (int)floor(-log(Precision)/log(2*D/Rmax)) ;
    if (NbrDir > 20){ if (Flag_FMM == 1) Msg(WARNING,"NbrDir = %d changed to 20 for memory requirements", NbrDir) ; NbrDir = 20 ;}
  }
 

  if (Flag_FMM == 2)  
    NbrDir = (NbrDir > Current.FMM.NbrDir) ? Current.FMM.NbrDir : NbrDir ; 
  else
    Current.FMM.NbrDir = (NbrDir>Current.FMM.NbrDir) ? NbrDir : Current.FMM.NbrDir ;
  
  GetDP_Return(NbrDir);
}



int  NeighbouringGroups( int FMMGroupEObs, int FMMGroupESrc ){

  int iFMMEqu, NbrFMMEqu, FMMObs, FMMSrc, NbrNG=0, i, *NG ;
  struct FMMmat  *FMMmat_P0 ;
  List_T *NG_L ;

  GetDP_Begin("NeigbouringGroups");

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0) ;
  
  iFMMEqu = 0 ; 
  while (iFMMEqu<NbrFMMEqu && ((FMMmat_P0+iFMMEqu)->Obs != Current.FMM.Obs || (FMMmat_P0+iFMMEqu)->Src != Current.FMM.Src)) 
    iFMMEqu++ ;

  if ( iFMMEqu == NbrFMMEqu )  Msg(GERROR, "Wrong Supports for FMM: Source %d Observation %d", Current.FMM.Src, Current.FMM.Obs) ; 

  FMMObs = (FMMmat_P0+iFMMEqu)->Obs ;
  FMMSrc = (FMMmat_P0+iFMMEqu)->Src ;

  if (FMMObs == FMMSrc &&  FMMGroupEObs == FMMGroupESrc) GetDP_Return(1) ;


  List_Read((FMMmat_P0+iFMMEqu)->NG_L, FMMGroupESrc, &NG_L);
  NbrNG = List_Nbr(NG_L) ;  

  NG = (int*)(NG_L->array) ;   
  
  for ( i = 0 ; i < NbrNG ; i++ )
    if (NG[i] == FMMGroupEObs) GetDP_Return(1) ; 
     
  GetDP_Return(0) ;
}



int  Get_NextElementSourceNeighbour( struct Element      * Element ) {

  struct Element * ElementSource ;

  GetDP_Begin("Get_NextElementSourceNeighbour");
  
  ElementSource = Element->ElementSource ;

  while (Get_NextElementSource(ElementSource))
    if (NeighbouringGroups(Element->FMMGroup, ElementSource->FMMGroup)) GetDP_Return(1) ;
  
  GetDP_Return(0) ;
}


int  Get_NextElementSourceInGroup(struct Element *Element) {
  struct Element * ElementSource ;

  GetDP_Begin("Get_NextElementSourceInGroup");

  ElementSource = Element->ElementSource ;
  while (Get_NextElementSource(ElementSource))
    if (ElementSource->FMMGroup == Element->FMMGroup) GetDP_Return(1) ;  
  GetDP_Return(0) ;
}



void  Geo_CreateFMMGroup( int InSupport, struct GeoData *GeoData_P, double k0 ){

  int NbrGroupsInSupport, iGroup, NbrFMMGroupsInRegion ;
  int NbrElms, NbrElmsInRegion, NbrElmsInSupport ;
  int j, k, l, m, FoundHome, iFMMGroup, iElm ;
  int ***NbrElmsDiv, ****ElmsDiv, Dimension ;
  double  Xmax, Xmin, Ymax, Ymin, Zmax, Zmin, l0 ;
  int NbrDivX = 1, NbrDivY = 1, NbrDivZ = 1 ;
  double  x1, x2, y1, y2, z1, z2, Centroid[3], tol ;
  struct Value DivXYZ ;
  struct Geo_Element *Elements_P0 ;
  struct FMMData FMMData_S ;
  struct FMMGroup FMMGroup_S ;
  List_T * FMMGroup_L ;  


  GetDP_Begin("Geo_CreateFMMGroup");

  /* Dimension = (int)GeoData_P->Dimension ; */
  Dimension = Current.FMM.Dimension ;    

  tol = 1e-6;

  if (!List_Nbr(Problem_S.FMMGroup)){
    if(k0 > 0){
      l0 = k0/TWO_PI ; /* Wavelength */
      Msg(INFO,"Helmholtz %dD, Wavelength = %g", Dimension, l0);
    }  
    else{
      Msg(INFO,"Laplace %dD", Dimension);
    }
  }

  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)) ; 
  
  NbrElmsInSupport = 0 ;  
  NbrGroupsInSupport = List_Nbr(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList ); 

  iFMMGroup = 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) ;
    NbrDivX = (int)DivXYZ.Val[0] ;
    NbrDivY = (int)DivXYZ.Val[1] ;
    NbrDivZ = (int)DivXYZ.Val[2] ;
  
    Msg(INFO,"Region = %d DivX = %d  DivY = %d  DivZ = %d", Current.Region, NbrDivX, NbrDivY, NbrDivZ) ; 
    
    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));
      }
    }

    for (iElm = 0 ; iElm < NbrElms ; iElm++){