fmm_operations.c

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

	  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->NumEqur, 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_AllLaplace3D(struct FMMmat *FMMmat_P, double *x, double *y ){
  /* Aggregation and Disaggregation matrices 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, l, m, j, k ;
  int  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_AllLaplace3D");

  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-1)*(Current.FMM.NbrDir+1) ;       
      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]-1)*(Nd_A[iG2]+1) ;
	  List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
	  TG1G2 = T[iG1][iG2] ;
	  
	  l = 0 ; m = 0 ; 
	  for (iDir = 0; iDir < NbrDir ; iDir++){
	    j = 0 ; k = 0 ;
	    TARe[iDir] = TAIm[iDir] = 0 ;	     
	    for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++){
	      iR = 2*((l+j)*(l+j+1)+m+k) ;
	      Treal = TG1G2[iR  ] ;
	      Timag = TG1G2[iR+1] ;	
	      
	      TARe[iDir] += Treal * AgJRe[iDir2] - Timag * AgJIm[iDir2] ;
	      TAIm[iDir] += Treal * AgJIm[iDir2] + Timag * AgJRe[iDir2] ;
	      
	      if (j==k) {j++;k=-j;} else k++;
	    }/* iDir2 NbrDir */
	    if (l==m) {l++;m=-l;} else m++;
	  }/* iDir NbrDir */

	  if ( (TypeTimeDerivative == NODT_) ||  (TypeTimeDerivative == DT_ && NbrHar == 1) ){
	    if(iHar == 0)
	      List_Read(FMMmat_P->NumEqur, NumFG[iG2], &NumEqu_L ) ; 
	    else
	      List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;  
	  }
	  
	  if (TypeTimeDerivative == DT_ && NbrHar == 2){
	    if (iHar == 0)
	      List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;
	    else
	      List_Read(FMMmat_P->NumEqur, 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)		
	      if ((TypeTimeDerivative == NODT_)|| (TypeTimeDerivative == DT_ && NbrHar == 1))
		y[NumEqu] +=  Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir] ;
	      else
		if (iHar ==0)		
		  y[NumEqu] += w * (Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir]) ;
		else
		  y[NumEqu] -= w * (Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir]) ;
		
	    
	  }/* iEqu NbrEqu */
	  
	}/* iG2 Destination */
      }/* if NbrFG */
    }/* iG1 Origin */    
  }/* iHar */
  
  GetDP_End ;
  
}



void FMMProd_AllLaplace3D_VECTOR(struct FMMmat *FMMmat_P, double *x, double *y ){
  /* Aggregation and Disaggregation matrices can be either SCALAR or VECTOR*/
  int  NbrGroupSrc, NbrDir ;
  int  NbrFG,  *NumFG, *Nd_A ; 
  int  iDof, NbrDof, NumDof, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, NumEquY, *NumEqu_A, *NumEquY_A ;
  int  iG1, iG2, iDir, iDir2, iR, iHar, NbrHar, l, m, j, k, iCom, NbrCom,Inc,iA,iTA ;
  int  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 ; 
  List_T *FG_L, *Nd_L, *NumDof_L, *NumEqu_L, *NumEquY_L ;


  GetDP_Begin("FMMProd_AllLaplace3D");


  T = FMMmat_P->T ;
  
  TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
  w = FMMmat_P->Pulsation ;

  NbrHar = Current.NbrHar ;
  NbrCom = FMMmat_P->NbrCom ;
  Inc = 2*NbrCom;

  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-1)*(Current.FMM.NbrDir+1) ;       
      for (iDir = 0; iDir < NbrCom* 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->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+=Inc){
	    for (iCom = 0 ; iCom < NbrCom ; iCom++){
	      AgJRe[iDir*NbrCom+iCom] += Ag_V[iR+iCom] * x[NumDof] ;
	      AgJIm[iDir*NbrCom+iCom] += Ag_V[iR+NbrCom+iCom] * x[NumDof] ;	  
	    }
	  }  
	} /*iDof NbrDof*/
      
	NumFG = (int*)(FG_L->array) ;
      
	for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
	  NbrDir = (Nd_A[iG2]-1)*(Nd_A[iG2]+1) ;
	  List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
	  TG1G2 = T[iG1][iG2] ;
		
	  l = 0 ; m = 0 ; 
	  for (iDir = 0; iDir < NbrDir ; iDir++){
	    j = 0 ; k = 0 ;
	    TARe[iDir] = TAIm[iDir] = 0 ;	     
	    for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++){
	      iR = 2*((l+j)*(l+j+1)+m+k) ;
	      Treal = TG1G2[iR  ] ;
	      Timag = TG1G2[iR+1] ;	
	      
	      for (iCom = 0 ; iCom < NbrCom ; iCom++){
		iTA = iDir*NbrCom+iCom ;
		iA = iDir2*NbrCom+iCom ;
		TARe[iTA] += Treal * AgJRe[iA] - Timag * AgJIm[iA] ;
		TAIm[iTA] += Treal * AgJIm[iA] + Timag * AgJRe[iA] ;
	      }
	      if (j==k) {j++;k=-j;} else k++;
	    }/* iDir2 NbrDir */
	    if (l==m) {l++;m=-l;} else m++;
	  }/* iDir NbrDir */

	  if(iHar == 0)
	    List_Read(FMMmat_P->NumEqu, NumFG[iG2], &NumEqu_L ) ; 
	  else
	    List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;  

	  if (TypeTimeDerivative == NODT_)
	    NumEquY_L = NumEqu_L ;
	  else
	    if (iHar == 0)
	      List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEquY_L ) ;
	    else
	      List_Read(FMMmat_P->NumEqu, NumFG[iG2], &NumEquY_L ) ;

	  NbrEqu = List_Nbr( NumEqu_L ) ;
	  NumEqu_A = (int*)(NumEqu_L->array) ;
	  NumEquY_A = (int*)(NumEquY_L->array) ;
      
	  for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	    NumEqu =  NumEqu_A[iEqu] - 1 ;
	    NumEquY =  NumEquY_A[iEqu] - 1 ;
	    Disag_V = Disag_M[iEqu] ;	
	    for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc)
	      for (iCom = 0 ; iCom < NbrCom ; iCom++){
		iTA = iDir*NbrCom+iCom ;
		
		if (TypeTimeDerivative == NODT_)
		  y[NumEqu] +=  Disag_V[iR+iCom] * TARe[iTA] - Disag_V[iR+NbrCom+iCom] * TAIm[iTA] ;
		else
		  if (iHar ==0)		
		    y[NumEquY] -= w * (Disag_V[iR+iCom] * TARe[iTA] - Disag_V[iR+NbrCom+iCom] * TAIm[iTA]) ;
		  else
		    y[NumEquY] += w * (Disag_V[iR+iCom] * TARe[iTA] - Disag_V[iR+NbrCom+iCom] * TAIm[iTA]) ;


	      }	     
	  }/* iEqu NbrEqu */
	  
	}/* iG2 Destination */
      }/* if NbrFG */
    }/* iG1 Origin */    
  }/* iHar */
  
  GetDP_End ;

}


void FMMProd_AllHelmholtz(struct FMMmat *FMMmat_P, double *x, double *y ){

  int  NbrGroupSrc, NbrDir; 
  int  NbrFG, *NumFG ;
  int  iDof, NbrDof, NumDof, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, NumEqui, *NumEqu_A, *NumEqui_A ;
  int  iG1, iG2, iDir, iHar,  NbrHar, iR, iCom, NbrCom,iA,iRD ;

  double  **Ag_M, *Ag_V, AgJRe[3*NBR_MAX_DIR], AgJIm[3*NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
  double TrAr, TrAi, TiAi, TiAr ;

  List_T *FG_L, *NumDof_L, *NumEqu_L,*NumEqui_L ;

  GetDP_Begin("FMMProd_AllHelmholtz");

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

  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++)
	for (iCom = 0 ; iCom < NbrCom ; iCom++)
	  AgJRe[iDir*NbrCom+iCom] = AgJIm[iDir*NbrCom+iCom] = 0 ;

      List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
      NbrFG = List_Nbr(FG_L) ;  

      if(NbrFG){      
	if(iHar == 0)
	  List_Read(FMMmat_P->NumDofr,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*NbrCom){
	    for (iCom = 0 ; iCom < NbrCom ; iCom++){
	      if(iHar == 0){
		AgJRe[iDir*NbrCom+iCom] += Ag_V[iR+iCom  ]      * x[NumDof];  
		AgJIm[iDir*NbrCom+iCom] += Ag_V[iR+NbrCom+iCom] * x[NumDof];  
	      }
	      else{
		AgJRe[iDir*NbrCom+iCom] -= Ag_V[iR+NbrCom+iCom] * x[NumDof];  
		AgJIm[iDir*NbrCom+iCom] += Ag_V[iR+iCom]   * x[NumDof]; 
	      }
	    }
	  }
	}

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

	  List_Read(FMMmat_P->NumEqur, NumFG[iG2], &NumEqu_L) ;
	  List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqui_L) ;

	  NbrEqu = List_Nbr( NumEqu_L ) ;
	  NumEqu_A = (int*)( NumEqu_L->array ) ;
	  NumEqui_A = (int*)( NumEqui_L->array ) ;

	  for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	    NumEqu =  NumEqu_A[iEqu] - 1 ;
	    NumEqui =  NumEqui_A[iEqu] - 1 ;
	    Disag_V = Disag_M[iEqu] ;

	    for (iDir = 0, iR = 0, iRD = 0 ; iDir < NbrDir ; iDir++, iR +=2, iRD+=2*NbrCom){
	      for (iCom = 0 ; iCom < NbrCom ; iCom++){
		iA = iDir* NbrCom +iCom ;
		TrAr = TG1G2[iR]   * AgJRe[iA] ;
		TrAi = TG1G2[iR]   * AgJIm[iA] ;
		TiAi = TG1G2[iR+1] * AgJIm[iA] ;
		TiAr = TG1G2[iR+1] * AgJRe[iA] ;
		y[NumEqu]  += Disag_V[iRD+iCom] * ( TrAr - TiAi ) - Disag_V[iRD+NbrCom+iCom] * ( TrAi + TiAr ) ;	
		y[NumEqui] += Disag_V[iRD+iCom] * ( TrAi + TiAr ) + Disag_V[iRD+NbrCom+iCom] * ( TrAr - TiAi ) ; 
	      }
	    }
	  }/* for iEqu NbrEqu*/
	}/* for iG2 Destination*/
      }/* if NbrFG */
    }/* for iG1 Origin */    
  }/* for iHar */

  GetDP_End ;  
}