fmm_groups.c

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

#define RCSID "$Id: FMM_Groups.c,v 1.17 2006/02/26 00:42:53 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 *
 * Please report all bugs and problems to <[e-mail]>.
 *
 * Contributor(s):
 *   Ruth Sabariego
 */

#include "GetDP.h"
#include "Numeric.h"
#include "Data_DefineE.h"
#include "CurrentData.h"
#include "Tools.h"
#include "Quadrature.h"
#include "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "Data_FMM.h"
#include "F_FMM.h"
#include "Truncation.h"


void Geo_InitFMMData( struct FMMData *FMMData_P ){
 
  GetDP_Begin("Geo_InitFMMData") ;
  
  FMMData_P->Group = -1 ;
  FMMData_P->Xgc = FMMData_P->Ygc = FMMData_P->Zgc = 0;
  FMMData_P->Rmax = -1. ;
  FMMData_P->Element = NULL ;  
  
  GetDP_End;
}


void InitFMMmatrix( int i_Src, int i_Obs, struct FMMmat *FMMmat_P ){

  GetDP_Begin("InitFMMmatrix") ;

  FMMmat_P->Src = i_Src ;
  FMMmat_P->Obs = i_Obs ;
  
  FMMmat_P->FunctionSpaceIndexDof = FMMmat_P->FunctionSpaceIndexEqu = -1 ;
  FMMmat_P->TypeTimeDerivative = NODT_ ; /* Default: No time derivative */
 
  FMMmat_P->Pulsation = 1. ;
  
  FMMmat_P->NbrCom = 0 ;

  FMMmat_P->NumDof = FMMmat_P->NumEqu = NULL ;
  FMMmat_P->NumDofr = FMMmat_P->NumEqur = NULL ; /* For Renumbering */
  FMMmat_P->NumDofi = FMMmat_P->NumEqui = NULL ; /* Necessary for renumbering when complex system*/
 
  FMMmat_P->NG_L = FMMmat_P->FG_L = NULL ;
  FMMmat_P->Nd_L = NULL ;

  FMMmat_P->A_L = FMMmat_P->D_L = NULL;
  FMMmat_P->T   = NULL ;
  FMMmat_P->FctProd = NULL ;

  GetDP_End;
}


void Init_CurrentFMMData ( double k0 ){

  int i, i_Phi, i_The, NbrDir, Val_ns=0, Dimension ;
  
  GetDP_Begin("Init_CurrentFMMData") ;


  Dimension = Current.FMM.Dimension ;

  if (k0 < 0){
    NbrDir = Current.FMM.NbrDir ;
  } 
 else {
    Val_ns = FMM_NbrSpectralDirections( k0 );
    NbrDir = (Dimension ==_2D) ? 2*Val_ns+1 : 2*Val_ns*Val_ns ;
    Current.FMM.NbrDir = NbrDir ; 
    Msg(INFO, "NbrDirections = %d Val_ns = %d", NbrDir, Val_ns) ;  
 }


  Current.FMM.N = (k0 < 0 && Dimension == _3D) ? (NbrDir-1)*(NbrDir+1): Current.FMM.NbrDir ; /* Laplace & GradLaplace 3D */

  if (Current.FMM.N > NBR_MAX_DIR || Current.FMM.NbrDir >  NBR_MAX_DIR )
    Msg(GERROR, "NbrDirections exceeds the NBR_MAX_DIR: %d", NBR_MAX_DIR ) ;

  if (k0 < 0) Msg(INFO, "NbrDirections = %d N = %d", NbrDir, Current.FMM.N ) ;  
 
  Current.FMM.Type    = SCALAR ; /* Default */
  Current.FMM.NbrCom  = 1 ;
  Current.FMM.Flag_GF = 0 ;
  Current.FMM.Xgc = Current.FMM.Ygc = Current.FMM.Zgc = 0. ;

  if (k0 < 0){
    Current.FMM.Phi    = NULL ;
    Current.FMM.Theta  = NULL ;
    Current.FMM.Weight = NULL ;
    Current.FMM.Kdir   = NULL ;
  }    
  else
    switch(Dimension){
    case _2D :
      Current.FMM.Phi = (double*)Malloc(NbrDir*sizeof(double));
      for (i_Phi = 0 ; i_Phi < NbrDir ; i_Phi++)  
	Current.FMM.Phi[i_Phi] = i_Phi* 2.*PI/NbrDir;    
      break;
      
    case _3D : 
      Current.FMM.Phi = (double*)Malloc(2*Val_ns*sizeof(double));
      Current.FMM.Theta  = (double*)Malloc(Val_ns*sizeof(double));
      Current.FMM.Weight = (double*)Malloc(NbrDir*sizeof(double));
      Current.FMM.Kdir = (double**)Malloc(NbrDir*sizeof(double));
      for (i = 0 ; i < NbrDir ; i++)
	Current.FMM.Kdir[i] = (double*)Malloc(3*sizeof(double));
      
      for (i_Phi = 0 ; i_Phi < 2*Val_ns ; i_Phi++)  
	Current.FMM.Phi[i_Phi] = i_Phi*PI/Val_ns ;   
      
      GaussLegendre(-1., 1. , Current.FMM.Theta, Current.FMM.Weight, Val_ns);
      for (i_The = 0 ; i_The < Val_ns ; i_The++)  
	Current.FMM.Theta[i_The] = acos(Current.FMM.Theta[i_The]);    
      
      for (i_Phi = 0 ; i_Phi < 2 * Val_ns ; i_Phi++){
	for (i_The=0 ; i_The < Val_ns ; i_The++){
	  Current.FMM.Kdir[i_Phi * Val_ns + i_The ][0] =
	    sin(Current.FMM.Theta[i_The])*cos(Current.FMM.Phi[i_Phi]);
	  Current.FMM.Kdir[i_Phi * Val_ns + i_The ][1] =
	    sin( Current.FMM.Theta[i_The])*sin(Current.FMM.Phi[i_Phi]);
	  Current.FMM.Kdir[i_Phi * Val_ns  + i_The ][2] =
	    cos(Current.FMM.Theta[i_The]);
	  Current.FMM.Weight[i_Phi * Val_ns  + i_The] =
	    Current.FMM.Weight[i_The];
	}
      }
      break;
    default :
      Msg(GERROR, "Wrong Dimension for GroupsFMM: %d", Dimension ) ;
    }
  
  GetDP_End;
  
}



void Get_RmaxGroups(int i_Obs, int i_Src){

  int Obs_OR_Src, iSup, NbrFMMSupports ; 
  int i, NbrFMMGroups ;
  struct FMMData      *FMMData_P0, *FMMData_P ;
  struct FMMGroup      FMMGroup_S ;

  GetDP_Begin("Get_RmaxGroups") ;
  
  NbrFMMSupports = (i_Obs == i_Src) ? 1 : 2 ;

  for(iSup = 0 ; iSup < NbrFMMSupports ; iSup ++ ){
    Obs_OR_Src = (iSup == 0) ? i_Obs :i_Src ;
    List_Read(Problem_S.FMMGroup, Obs_OR_Src, &FMMGroup_S) ;

    NbrFMMGroups = List_Nbr(FMMGroup_S.List) ;
    FMMData_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;
 
    for ( i = 0 ; i < NbrFMMGroups; i++){
      FMMData_P = FMMData_P0 + i;
      if(FMMData_P->Rmax == -1)
   	FMMData_P->Rmax = Geo_GetRmaxInFMMGroup(FMMData_P) ;   
    }/* NbrFMMGroups */
  }/*NbrFMMSupports */

  GetDP_End ;
}


double Geo_GetRmaxInFMMGroup( struct FMMData *FMMData_P){
 
  int iElm, i, NbrElmsGr, NumElmi ;
  double R , Rmax =0., Xgc, Ygc, Zgc ;
  double x1[NBR_MAX_NODES_IN_ELEMENT], y1[NBR_MAX_NODES_IN_ELEMENT], z1[NBR_MAX_NODES_IN_ELEMENT];
  struct Geo_Element *Element ;  

  GetDP_Begin("Geo_GetRmaxInFMMGroup");

  Xgc = FMMData_P->Xgc ;
  Ygc = FMMData_P->Ygc ;
  Zgc = FMMData_P->Zgc ;
  NbrElmsGr = List_Nbr(FMMData_P->Element) ; 
  
  for (iElm = 0 ; iElm < NbrElmsGr ; iElm++){
    List_Read(FMMData_P->Element, iElm, &NumElmi);
    Element =  Geo_GetGeoElementOfNum(NumElmi) ; 
    Geo_GetNodesCoordinates(Element->NbrNodes, Element->NumNodes,
			    x1, y1, z1);

    for (i = 0 ; i < Element->NbrNodes ; i++){    
      R = sqrt(SQU(x1[i]-Xgc)+ SQU(y1[i]-Ygc)+ SQU(z1[i]-Zgc)) ;  
      Rmax = ( R > Rmax ) ? R : Rmax ;
    }
  }

  GetDP_Return(Rmax) ;
}




void Get_GroupNeighbours( int i_FMMEqu ){

  int i_Obs, i_Src, Ndir ;
  int NbrGroupObs, NbrGroupSrc, i, j, Flag_Far = 0 ;
  double d, XSrc, YSrc, ZSrc, Dfar ;
  double **Mat =NULL ;
  struct FMMData      *FMMDataGObs_P0, *FMMDataGSrc_P0 ;
  struct FMMGroup      FMMGroup_S ;
  struct FMMmat       *FMMmat_P, *FMMmat_P0 ;
  List_T *NGi, *FGi, *Ndi ;

  GetDP_Begin("Geo_GroupNeighbours");


  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0) ;
  FMMmat_P  = FMMmat_P0 + i_FMMEqu ;
  i_Obs = FMMmat_P->Obs ;
  i_Src = FMMmat_P->Src ;

  Get_RmaxGroups(i_Obs, i_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) ;
    }

    FMMmat_P->A_L = List_Create(NbrGroupSrc, 1, sizeof(double**)) ;
    for(i = 0 ; i < NbrGroupSrc ; i++) List_Add(FMMmat_P->A_L, &Mat) ;

    FMMmat_P->D_L = List_Create(NbrGroupObs, 1, sizeof(double**)) ;
    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) ;
    }

    FMMmat_P->NG_L = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
    FMMmat_P->FG_L = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
    FMMmat_P->Nd_L = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;

    FMMmat_P->A_L = List_Create(NbrGroupSrc, 1, sizeof(double**)) ;
    for(i = 0 ; i < NbrGroupSrc ; i++) List_Add(FMMmat_P->A_L, &Mat) ;

    FMMmat_P->D_L = List_Create(NbrGroupObs, 1, sizeof(double**)) ;
    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)) ;

      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)) ;
       
	/* Dfar = Current.FMM.far * (FMMDataGSrc_P0+i)->Rmax *2 ; */
	Dfar = Current.FMM.far ;
	if(d < Dfar)
	  List_Add(NGi, &j) ;
	else{ 
	  List_Add(FGi, &j) ;

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

	  if (!Flag_Far) Flag_Far = 1;
	} 
      }/* NbrGroupObs */
      List_Add( FMMmat_P->NG_L, &NGi );
      List_Add( FMMmat_P->FG_L, &FGi );
      List_Add( FMMmat_P->Nd_L, &Ndi );
    }/* NbrGroupSrc */
  } 

  GetDP_End ;
}


void  Geo_SetCentroidCoordinates(int Num, double C[3]){
 
  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_SetCentroidCoordinates");
  
  C[0] = C[1] = C[2] = 0.;

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

  for (j = 0;  j < Element->NbrNodes ; j++){
    C[0] += x[j] ;
    C[1] += y[j] ;
    C[2] += z[j] ;
  }
  
  C[0] /= Element->NbrNodes ;
  C[1] /= Element->NbrNodes ;
  C[2] /= Element->NbrNodes ;
  
  GetDP_End ;
}


void ReSet_FMMGroupCenters( ){

  int i, NbrGroups, iG, NbrElms, iElm, NumElm ;
  double Centroid[3];
  struct FMMGroup     FMMGroup_S ;
  struct FMMData     *FMMDataG_P0, *FMMDataG_P ;


  GetDP_Begin("ReSet_FMMGroupCenters");
  
  for( i = 0 ; i < List_Nbr(Problem_S.FMMGroup); i++){
    List_Read(Problem_S.FMMGroup, i, &FMMGroup_S) ;
    NbrGroups = List_Nbr(FMMGroup_S.List) ;
    FMMDataG_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;
    
    for( iG = 0 ; iG < NbrGroups; iG++){
      FMMDataG_P =  FMMDataG_P0 + iG ;
      FMMDataG_P->Xgc =  FMMDataG_P->Ygc = FMMDataG_P->Zgc = 0. ;      
      NbrElms = List_Nbr(FMMDataG_P->Element) ;

      for( iElm = 0 ; iElm < NbrElms ; iElm++){
	List_Read( FMMDataG_P->Element, iElm, &NumElm) ;
	Geo_SetCentroidCoordinates(NumElm, Centroid);
	FMMDataG_P->Xgc += Centroid[0] ;  
	FMMDataG_P->Ygc += Centroid[1] ;
	FMMDataG_P->Zgc += Centroid[2] ;
      }/* iElm NbrElms */

      FMMDataG_P->Xgc /= NbrElms ;  
      FMMDataG_P->Ygc /= NbrElms ;
      FMMDataG_P->Zgc /= NbrElms ;
    }/* iG FMMGroups */

  }/* i  Getdp Groups affected by FMM */

  GetDP_End;

}


void ReGenerate_FMMGroupNeighbours(int Flag_FMMDA){

  int i_Obs, i_Src, Ndir, FMMSrc,FMMObs ;
  int NbrGroupObs, NbrGroupSrc, i, j, NbrFMMEqu, i_FMMEqu, i_Group  ;
  int  Nbr_DofList, Nbr_EquList ;
  double d, XSrc, YSrc, ZSrc, Dfar ;
  double **Mat =NULL, **Aggreg_M, **Disagg_M ;
  struct FMMData      *FMMDataGObs_P0, *FMMDataGSrc_P0 ;
  struct FMMGroup      FMMGroup_S ;
  struct FMMmat       *FMMmat_P, *FMMmat_P0 ;
  List_T *NGi, *FGi, *Ndi, * NumDof_L, * NumEqu_L ;


  GetDP_Begin("ReGenerate_FMMGroupNeighbours");

  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0) ;
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix);
  Dfar = Current.FMM.far ;

  /* With movement neighbouring and far groups have to be recalculated */

  for (i_FMMEqu = 0 ; i_FMMEqu < NbrFMMEqu ; i_FMMEqu++){
    FMMmat_P  = FMMmat_P0 + i_FMMEqu ;
    List_Reset(FMMmat_P->NG_L);
    List_Reset(FMMmat_P->FG_L);
    List_Reset(FMMmat_P->Nd_L);
    
    if(Flag_FMMDA){
      FMMSrc = (FMMmat_P0+i_FMMEqu)->Src ;
      List_Read(Problem_S.FMMGroup, FMMSrc, &FMMGroup_S ) ;   
      NbrGroupSrc  = List_Nbr(FMMGroup_S.List ) ;
      for ( i_Group = 0 ; i_Group < NbrGroupSrc ; i_Group++ ) {
	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) ;
	for (j = 0 ; j < Nbr_DofList ; j++)
	  Free(Aggreg_M[j]) ;
	Free(Aggreg_M);
      }

      FMMObs = (FMMmat_P0+i_FMMEqu)->Obs ;
      List_Read(Problem_S.FMMGroup, FMMObs, &FMMGroup_S ) ;   
      NbrGroupObs  = List_Nbr(FMMGroup_S.List ) ;

      for ( i_Group = 0 ; i_Group < NbrGroupObs ; i_Group++ ) {
	List_Read(FMMmat_P->D_L, i_Group, &Disagg_M) ;
	List_Read(FMMmat_P->NumEqu, i_Group, &NumEqu_L) ;
	Nbr_EquList = List_Nbr(NumEqu_L) ;
	for (j = 0 ; j < Nbr_EquList ; j++)
	  Free(Disagg_M[j]);
	Free(Disagg_M);