fmm_operations.c

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

  List_T *NumDof_L, *NumEqu_L ;


  GetDP_Begin("FMM_InverseRenumbering");

  /* This function does not work properly! To revise... */

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;
    NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
  
    for (iG = 0 ; iG < NbrGroupSrc ; iG++){
      List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
      NbrDof = List_Nbr( NumDof_L);   
      NumDof_A = (int*)(NumDof_L->array) ;
      
      for (iDof = 0 ; iDof < NbrDof ; iDof++){
	NumDofr  = NumDof_A[iDof] ;
	NumDof = rpermr[NumDofr]-1 ;
	List_Write(NumDof_L, iDof, &NumDof) ;
      }
      List_Write(FMMmat_P->NumDof, iG, &NumDof_L ) ;
    }
 
    if(FMMmat_P->NumEqu !=  FMMmat_P->NumDof){
      NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
 
      for (iG = 0 ; iG < NbrGroupObs ; iG++){
	List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
	NbrEqu = List_Nbr( NumEqu_L);   
	NumEqu_A = (int*)(NumEqu_L->array) ;

	for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	  NumEqur = NumEqu_A[iEqu] ;
	  NumEqu  = rpermr[NumEqur]-1 ;
	  List_Write(NumEqu_L, iEqu, &NumEqu ) ;
	}
	List_Write(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
      }
    }
  }

 
  GetDP_End;
}






void FMM_Scaling(double *rowscal, double *colscal, int N){

  int  NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs;
  int  iR, iDir, NbrDir, iCom, NbrCom, Inc, TypeTimeDerivative, i ,Flag_scaled = 0 ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  int  iDof, NbrDof, NumDof, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, *NumEqu_A ;
  double **Ag_M, **D_M, *Ag_V, *D_V, *x_prev = NULL ;
  List_T *NumDof_L, *NumEqu_L ;

  /* Aggregation and Disaggregation matrices can be either SCALAR or VECTOR */

  GetDP_Begin("FMM_Scaling");
  
  
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;

  NbrDir = Current.FMM.N ; 


  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;
    
    TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
   
    if ( TypeTimeDerivative == DT_ && Current.NbrHar == 1 ){
#if defined(HAVE_SPARSKIT)
      x_prev = ((Current.DofData->CurrentSolution-1)->x).V ;
#else
      Msg(GERROR, "FMM_Scaling only implemented for Sparskit");
#endif
      if (!Flag_scaled){
	for (i = 0 ; i < N ; i++)
	  x_prev[i] /= colscal[i] ;
	Flag_scaled = 1;
      }
    }
    
    NbrCom = FMMmat_P->NbrCom ;
    Inc = 2 * NbrCom ;
    NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
    
    for (iG = 0 ; iG < NbrGroupSrc ; iG++){
      List_Read(FMMmat_P->A_L, iG, &Ag_M ) ;
      List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
      NbrDof = List_Nbr( NumDof_L);   
      NumDof_A = (int*)(NumDof_L->array) ;      

      for (iDof = 0 ; iDof < NbrDof ; iDof++){
	NumDof  = NumDof_A[iDof] - 1 ;
	Ag_V = Ag_M[iDof] ;
	
	for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){
	  for (iCom = 0 ; iCom < NbrCom ; iCom++){
	    Ag_V[iR+iCom]   *= colscal[NumDof] ;
	    Ag_V[iR+NbrCom+iCom] *= colscal[NumDof] ;
	  }
	}	
      }
    }/* iG Source */
 
    NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
  
    for (iG = 0 ; iG < NbrGroupObs ; iG++){
      List_Read(FMMmat_P->D_L, iG, &D_M ) ;
      List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
      NbrEqu = List_Nbr( NumEqu_L);   
      NumEqu_A = (int*)(NumEqu_L->array) ;

      for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	NumEqu  = NumEqu_A[iEqu] - 1 ;
	D_V = D_M[iEqu] ;
	
	for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){	  
	  for (iCom = 0 ; iCom < NbrCom ; iCom++){
	    D_V[iR+iCom]   *= rowscal[NumEqu] ;
	    D_V[iR+NbrCom+iCom] *= rowscal[NumEqu] ;
	  }
	}	
      }
    }/* iG Observation */

  }/* iFMMEqu */
 
  GetDP_End;
}



void FMM_UnScaling(double *rowscal, double *colscal){

  int  NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs, iR, iDir, NbrDir, iCom, NbrCom, Inc ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  int  iDof, NbrDof, NumDof, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, *NumEqu_A ;
  double **Ag_M, **D_M, *Ag_V, *D_V ;
  List_T *NumDof_L, *NumEqu_L ;

  /* Aggregation and Disaggregation matrices can be either SCALAR or VECTOR*/

  GetDP_Begin("FMM_UnScaling");
  
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  NbrDir = Current.FMM.N ; 
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;

    NbrCom = FMMmat_P->NbrCom ;
    Inc = 2 * NbrCom ;
    NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
  
    for (iG = 0 ; iG < NbrGroupSrc ; iG++){
      List_Read(FMMmat_P->A_L, iG, &Ag_M ) ;
      List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
      NbrDof = List_Nbr( NumDof_L);   
      NumDof_A = (int*)(NumDof_L->array) ;
      
      for (iDof = 0 ; iDof < NbrDof ; iDof++){
	NumDof  = NumDof_A[iDof] - 1 ;
	Ag_V = Ag_M[iDof] ;

	for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){	 
	  for (iCom = 0 ; iCom < NbrCom ; iCom++){
	    Ag_V[iR+iCom]   /= colscal[NumDof] ;
	    Ag_V[iR+NbrCom+iCom] /= colscal[NumDof] ;
	  }
	}	
      }
    }/* iG Source */
 
    NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
  
    for (iG = 0 ; iG < NbrGroupObs ; iG++){
      List_Read(FMMmat_P->D_L, iG, &D_M ) ;
      List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
      NbrEqu = List_Nbr( NumEqu_L);   
      NumEqu_A = (int*)(NumEqu_L->array) ;

      for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	NumEqu  = NumEqu_A[iEqu] - 1 ;
	D_V = D_M[iEqu] ;
	for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){	
	  for (iCom = 0 ; iCom < NbrCom ; iCom++){
	    D_V[iR+iCom]   /= rowscal[NumEqu] ;
	    D_V[iR+NbrCom+iCom] /= rowscal[NumEqu] ;
	  }
	}	
      }
    }/* iG Observation */
    
  }/* iFMMEqu */
 
  GetDP_End;
}




void FMM_DTAMatrix(int N, double ***DTA){ 

  int  i, NbrFMMEqu, iFMMEqu ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  double **M, *x, *y ;

  GetDP_Begin("FMM_DTAMatrix");

  M =(double**)Malloc(N *sizeof(double*));
  for(i = 0 ; i < N ; i++)
    M[i] =(double*)Calloc(N,sizeof(double));

  x = (double*)Calloc(N, sizeof(double));

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;
    for (i = 0 ; i < N ; i++){
      x[i] = 1. ;
      y = M[i] ; /* row */
      ((void (*)(struct FMMmat*,double*,double*))FMMmat_P->FctProd)(FMMmat_P, x, y) ;
      x[i] = 0. ;
    }
  }
  
  *DTA = M ;
  GetDP_End;
}



/* FMMProd_Laplace2D and FMMProd_Laplace3D admit only SCALAR Aggregation and Disaggregation,
   the VECTOR case is not yet implemented and will be considered in separate routines */ 

void FMMProd_AllLaplace2DNonAdaptive(struct FMMmat *FMMmat_P, double *x, double *y ){
  /* WARNING: Aggregation and Disaggregation matrices must be SCALAR */
  int  NbrGroupSrc, NbrDir ;
  int  NbrFG, *NumFG ;
  int  iDof, NbrDof, NumDof, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, *NumEqu_A ;
  int  iG1, iG2, iDir, iDir2, iR, iHar, NbrHar ;
  double  **Ag_M, *Ag_V, AgJRe[NBR_MAX_DIR], AgJIm[NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
  double Treal, Timag, TARe[NBR_MAX_DIR], TAIm[NBR_MAX_DIR] ; 

  List_T *FG_L, *NumDof_L, *NumEqu_L ;


  GetDP_Begin("FMMProd_AllLaplace2DNonAdaptive");

  T = FMMmat_P->T ;
  
  NbrHar = Current.NbrHar ;
  NbrDir = Current.FMM.NbrDir ;

  NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
    
  for (iHar = 0 ; iHar < NbrHar ; iHar++){
    for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
      List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
      for (iDir = 0; iDir < NbrDir ; iDir++)
	AgJRe[iDir] = AgJIm[iDir] = 0.;
      
      List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
      NbrFG = List_Nbr(FG_L) ;  
      if(NbrFG){
	if (iHar == 0)
	  List_Read(FMMmat_P->NumDof, iG1, &NumDof_L ) ;
	else 
	  List_Read(FMMmat_P->NumDofi, iG1, &NumDof_L ) ;

	NbrDof = List_Nbr(NumDof_L);
	NumDof_A = (int*)(NumDof_L->array) ;
	
	for (iDof = 0 ; iDof < NbrDof ; iDof++){
	  NumDof = NumDof_A[iDof] - 1 ;
	  Ag_V = Ag_M[iDof] ;
	  for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){	  
	    AgJRe[iDir] += Ag_V[iR  ] * x[NumDof] ;
	    AgJIm[iDir] += Ag_V[iR+1] * x[NumDof] ;
	  }
	}  /* iDof NbrDof */
	
	NumFG = (int*)(FG_L->array) ;     
	for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
	  List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
	  TG1G2 = T[iG1][iG2] ;
	  
	  for (iDir = 0, iR = 0 ; iDir < NbrDir ; iDir++){ /* Translation is not diagonal! */
	    TARe[iDir] = TAIm[iDir] = 0. ;
	    for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++, iR+=2){
	      Treal = TG1G2[iR  ] ;
	      Timag = TG1G2[iR+1] ;
	      TARe[iDir] += Treal * AgJRe[iDir2] - Timag * AgJIm[iDir2] ;
	      TAIm[iDir] += Treal * AgJIm[iDir2] + Timag * AgJRe[iDir2] ;        
	    }/* iDir2 */
	  }/* iDir1 */
	  
	  if(iHar == 0)
	    List_Read(FMMmat_P->NumEqu, NumFG[iG2], &NumEqu_L ) ;
	  else
	    List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;

	  NbrEqu = List_Nbr( NumEqu_L ) ;
	  NumEqu_A = (int*)( NumEqu_L->array) ;
	  for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	    NumEqu =  NumEqu_A[iEqu] - 1 ;
	    Disag_V = Disag_M[iEqu] ;
	    for (iDir = 0, iR = 0 ; iDir < NbrDir ; iDir++, iR+=2)
	      y[NumEqu] += Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir] ;	     			   
	  }/* iEqu NbrEqu*/
	  
	}/* iG2 Destination*/
      }/* if NbrFG */
    }/* iG1 Origin */
  }/* iHar */
  
  GetDP_End ;
  
} 


void FMMProd_AllLaplace2D(struct FMMmat *FMMmat_P, double *x, double *y ){
  /* WARNING: Aggregation and Disaggregation matrices must be SCALAR */
  int  NbrGroupSrc, NbrDir ;
  int  NbrFG, *NumFG, *Nd_A ;
  int  iDof, NbrDof, NumDof, NumDofprev, *NumDof_A, *NumDofprev_A = NULL;
  int  iEqu, NbrEqu, NumEqu, *NumEqu_A ;
  int  iG1, iG2, iDir, iDir2, iR, iHar, NbrHar,TypeTimeDerivative ;
  double  **Ag_M, *Ag_V, AgJRe[NBR_MAX_DIR], AgJIm[NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
  double Treal, Timag, TARe[NBR_MAX_DIR], TAIm[NBR_MAX_DIR], w, *x_prev = NULL; 

  List_T *FG_L, *Nd_L, *NumDof_L, *NumDofprev_L, *NumEqu_L ;


  GetDP_Begin("FMMProd_AllLaplace2D");

  T = FMMmat_P->T ;
  
  NbrHar = Current.NbrHar ;
  
  TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
  w = FMMmat_P->Pulsation ;

  if (TypeTimeDerivative == DT_ && NbrHar == 1 )
#if defined(HAVE_SPARSKIT)
    x_prev = ((Current.DofData->CurrentSolution-1)->x).V ;
#else
    Msg(GERROR, "FMM_Scaling only implemented for Sparskit");
#endif


  NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
    
  for (iHar = 0 ; iHar < NbrHar ; iHar++){
    for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
      List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
      NbrDir = Current.FMM.NbrDir ;  
      for (iDir = 0; iDir < NbrDir ; iDir++)
	AgJRe[iDir] = AgJIm[iDir] = 0.;
      
      List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
      NbrFG = List_Nbr(FG_L) ;  
      if(NbrFG){
	List_Read(FMMmat_P->Nd_L, iG1, &Nd_L) ;
	Nd_A =  (int*)(Nd_L->array) ;
	
	if (iHar == 0)
	  List_Read(FMMmat_P->NumDofr, iG1, &NumDof_L ) ;
	else 
	  List_Read(FMMmat_P->NumDofi, iG1, &NumDof_L ) ;

	if(TypeTimeDerivative == DT_ && NbrHar ==1 ){
	  List_Read(FMMmat_P->NumDof, iG1, &NumDofprev_L ) ;
	  NumDofprev_A = (int*)(NumDofprev_L->array) ;
	}
	
	NbrDof = List_Nbr(NumDof_L);
	NumDof_A = (int*)(NumDof_L->array) ;
	
	for (iDof = 0 ; iDof < NbrDof ; iDof++){
	  NumDof = NumDof_A[iDof] - 1 ;
	  Ag_V = Ag_M[iDof] ;
	  
	  if(TypeTimeDerivative == NODT_ || ( TypeTimeDerivative == DT_ && NbrHar == 2) ){
	    for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){	  
	      AgJRe[iDir] += Ag_V[iR  ] * x[NumDof] ;
	      AgJIm[iDir] += Ag_V[iR+1] * x[NumDof] ;
	    }
	  }

	  if (TypeTimeDerivative == DT_ && NbrHar ==1 ){
	    NumDofprev = NumDofprev_A[iDof] - 1 ;
	    for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){	  
	      AgJRe[iDir] += Ag_V[iR  ] * ( x[NumDof] - x_prev[NumDofprev] ) ;
	      AgJIm[iDir] += Ag_V[iR+1] * ( x[NumDof] - x_prev[NumDofprev] ) ;
	    }
	  }

	 

	}  /* iDof NbrDof */
	
	NumFG = (int*)(FG_L->array) ;     
	for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
	  NbrDir = Nd_A[iG2] ;
	  List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;