print_fmm.c

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

	  fprintf(file,"%d(%d) ", FG[j], NbrDir);	

      if(NbrFG && Current.FMM.Precision != 0) NbrDirAV /=NbrFG ;
      fprintf(file, "\n\n");    
    }
    
    TotalMem += StorageDir + StorageFMM ;

  } /* iFMMEqu */
  
  fprintf(file,"\nAVERAGE Truncantion factor = %d\n", NbrDirAV) ;
  fprintf(file,"MAXIMUM Truncation factor = %d\n", NbrDirMAX);
  fprintf(file,"\nTotal Number of Dof = %d\n", Current.DofData->NbrAnyDof);
  fprintf(file,"TotalMemMoM %.5g Kb\n", Current.DofData->NbrAnyDof*Current.DofData->NbrAnyDof*32./1024.);
  fprintf(file, "\nTotalMemFMM %.5g Kb   StorageDir %.5g Kb   StorageFMM   %.5g Kb\n", TotalMem, StorageDir, StorageFMM);

  fprintf(file, "\nNbrMultLaplace2D %d\n", NbrMultLaplace2D) ;
  fprintf(file, "NbrMultLaplace3D %d\n", NbrMultLaplace3D) ;

  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu++){
    List_Read(Current.DofData->FMM_Matrix, iFMMEqu, &FMMmat_S ) ;
    List_Read(Problem_S.FMMGroup, FMMmat_S.Src, &FMMGroupSrc_S) ;
    FMMDataSrc_P0 = FMMDataObs_P0 =  (struct FMMData*)List_Pointer(FMMGroupSrc_S.List, 0) ; 
    NbrGroupsSrc = NbrGroupsObs = List_Nbr(FMMGroupSrc_S.List) ;

    if ( FMMmat_S.Src == FMMmat_S.Obs ){
      fprintf(file, "$EQUTERM (%d) SOURCE == OBSERVATION (%d)\n", iFMMEqu, FMMmat_S.Src) ;
      fprintf(file, "FunctionSpaceIndexDof %d FunctionSpaceIndexEqu %d\n",
	      FMMmat_S.FunctionSpaceIndexDof, FMMmat_S.FunctionSpaceIndexEqu) ;
    
      for (i = 0 ; i < NbrGroupsSrc ; i++){
	NbrElmsGroup =  List_Nbr((FMMDataSrc_P0+i)->Element) ;     
	fprintf(file, "$GROUP %d (Equ %d, Support %d)\n", i, iFMMEqu, FMMmat_S.Src) ;
	fprintf(file, "%f\n", (FMMDataSrc_P0+i)->Rmax); /* Group radius */
	fprintf(file, "%f %f %f\n", (FMMDataSrc_P0+i)->Xgc, (FMMDataSrc_P0+i)->Ygc, (FMMDataSrc_P0+i)->Zgc); /* Group Center */
      
	NbrElmsGroup =  List_Nbr((FMMDataSrc_P0+i)->Element) ;
	fprintf(file, "%d $EL ", NbrElmsGroup);
	for (j = 0 ; j < NbrElmsGroup ; j++){
	  List_Read((FMMDataSrc_P0+i)->Element, j, &NumElm);
	  fprintf(file, "%d ", NumElm);
	}
	fprintf(file, "\n");
	
	List_Read(FMMmat_S.NumDof, i, &NumDof_L) ;
	NbrDof = List_Nbr(NumDof_L) ;
	fprintf(file, "%d $Dof(%d)", NbrDof, FMMmat_S.FunctionSpaceIndexDof);
	NumDof = (int*)(NumDof_L->array) ;
	for (j = 0 ; j < NbrDof ; j++)
	  fprintf(file, "%d ", NumDof[j]);
	fprintf(file, "\n");

	List_Read(FMMmat_S.NumEqu, i, &NumEqu_L) ;
	NbrEqu = List_Nbr(NumEqu_L) ;
	fprintf(file, "%d $Equ(%d)", NbrEqu, FMMmat_S.FunctionSpaceIndexEqu);
	NumEqu = (int*)(NumEqu_L->array) ;
	for (j = 0 ; j < NbrEqu ; j++)
	  fprintf(file, "%d ", NumEqu[j]);
	fprintf(file, "\n");
	
	fprintf(file, "$END_GROUP %d \n", i);
      }
    }   
    else{ 
      List_Read(Problem_S.FMMGroup, FMMmat_S.Obs, &FMMGroupObs_S) ;
      FMMDataObs_P0 = (struct FMMData*)List_Pointer(FMMGroupObs_S.List, 0) ; 
      NbrGroupsObs = List_Nbr(FMMGroupObs_S.List) ;

      fprintf(file, "$EQUTERM (%d) SOURCE (%d) != OBSERVATION (%d)\n", iFMMEqu, FMMmat_S.Src, FMMmat_S.Obs ) ;
      fprintf(file, "FunctionSpaceIndexDof %d FunctionSpaceIndexEqu %d\n",
	      FMMmat_S.FunctionSpaceIndexDof, FMMmat_S.FunctionSpaceIndexEqu) ;
      fprintf(file, "$SOURCE GROUPS on %d (Equ %d) \n", FMMmat_S.Src, iFMMEqu) ;
      for (i = 0 ; i < NbrGroupsSrc ; i++){
	NbrElmsGroup =  List_Nbr((FMMDataSrc_P0+i)->Element) ;     
	fprintf(file, "$GROUP SRC %d (Equ %d)\n", i, iFMMEqu);
	fprintf(file, "%f\n", (FMMDataSrc_P0+i)->Rmax); /* Group radius */
	fprintf(file, "%f %f %f\n", (FMMDataSrc_P0+i)->Xgc, (FMMDataSrc_P0+i)->Ygc, (FMMDataSrc_P0+i)->Zgc); /* Group Center */
      
	NbrElmsGroup =  List_Nbr((FMMDataSrc_P0+i)->Element) ;
	fprintf(file, "%d $EL ", NbrElmsGroup);
	for (j = 0 ; j < NbrElmsGroup ; j++){
	  List_Read((FMMDataSrc_P0+i)->Element, j, &NumElm);
	  fprintf(file, "%d ", NumElm);
	}
	fprintf(file, "\n");
	List_Read(FMMmat_S.NumDof, i, &NumDof_L) ;
	NbrDof = List_Nbr(NumDof_L) ;
	fprintf(file, "%d $Dof(%d)", NbrDof, FMMmat_S.FunctionSpaceIndexDof);

	NumDof = (int*)(NumDof_L->array) ;
	for (j = 0 ; j < NbrDof ; j++)
	  fprintf(file, "%d ", NumDof[j]);
	fprintf(file, "\n");
	fprintf(file, "$END_GROUP SRC %d \n", i);
      }
      fprintf(file, "$END_SOURCE GROUPS on %d (Equ %d) \n", FMMmat_S.Src, iFMMEqu) ;

      fprintf(file, "$OBSERVATION GROUPS on %d (Equ %d) \n", FMMmat_S.Obs, iFMMEqu) ;
      for (i = 0 ; i < NbrGroupsObs ; i++){
	fprintf(file, "$GROUP OBS %d (Equ %d)\n", i, iFMMEqu);
	fprintf(file, "%f\n", (FMMDataObs_P0+i)->Rmax); /* Group radius */
	fprintf(file, "%f %f %f\n", (FMMDataObs_P0+i)->Xgc, (FMMDataObs_P0+i)->Ygc, (FMMDataObs_P0+i)->Zgc); /* Group Center */
	NbrElmsGroup =  List_Nbr((FMMDataObs_P0+i)->Element) ;
	fprintf(file, "%d $EL ", NbrElmsGroup);
	for (j = 0 ; j < NbrElmsGroup ; j++){
	  List_Read((FMMDataObs_P0+i)->Element, j, &NumElm);
	  fprintf(file, "%d ", NumElm);
	}
	fprintf(file, "\n");

	List_Read(FMMmat_S.NumEqu, i, &NumEqu_L) ;
	NbrEqu = List_Nbr(NumEqu_L) ;
	fprintf(file, "%d $Equ(%d)", NbrEqu, FMMmat_S.FunctionSpaceIndexEqu);	
	NumEqu = (int*)(NumEqu_L->array) ;
	for (j = 0 ; j < NbrEqu ; j++)
	  fprintf(file, "%d ", NumEqu[j]);
	fprintf(file, "\n");
	fprintf(file, "$END_GROUP OBS%d \n", i);
      }
      fprintf(file, "$END_OBSERVATION GROUPS on %d (Equ %d) \n", FMMmat_S.Obs, iFMMEqu) ;
    }
  }
    
  fclose(file);
  
  GetDP_End;
}


void Print_FMMGroupNeighbours( char* FileName ){
  
  int i, j, *NG, *FG, *Nd, NbrGroupsSrc ;
  int NbrNG, NbrFG ;
  int NbrFMMEqu, iFMMEqu ;
  struct FMMGroup FMMGroupSrc_S ;
  struct FMMmat FMMmat_S ;
  List_T *G_L, *Nd_L ;
 

  FILE * file ;
  

  GetDP_Begin("Print_FMMGroupNeighbours") ;
  
  Msg(INFO, "Writing FMM Group Neighbours : '%s' ", FileName );

  file = fopen(FileName, "w");
  fprintf(file, "$FMMGroupsNEIGHBOURS\n");
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  fprintf(file, "$NbrFMMEqu %d\n", NbrFMMEqu) ;

  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu++){
    List_Read(Current.DofData->FMM_Matrix, iFMMEqu, &FMMmat_S ) ;
    fprintf(file, "EQUTERM %d Src %d Obs %d FS_Dof %d FS_Equ %d\n",
	    iFMMEqu, FMMmat_S.Src, FMMmat_S.Obs, FMMmat_S.FunctionSpaceIndexDof, FMMmat_S.FunctionSpaceIndexEqu) ;

    List_Read(Problem_S.FMMGroup, FMMmat_S.Src, &FMMGroupSrc_S) ;    
    NbrGroupsSrc = List_Nbr(FMMGroupSrc_S.List) ;

    for ( i = 0 ; i < NbrGroupsSrc ; i++){
      fprintf(file, "FMMGroup %d (Equ %d on %d)", i, iFMMEqu, FMMmat_S.Src) ;
      List_Read(FMMmat_S.NG_L, i, &G_L) ;
      NbrNG = List_Nbr(G_L) ;
      NG = (int*)(G_L->array) ;
      fprintf(file, "\nNEIGHBOURS (on %d) ", FMMmat_S.Obs) ;
     
      for (j = 0 ; j < NbrNG ; j++ )
	fprintf(file, "%d ", NG[j]);

      List_Read(FMMmat_S.FG_L, i, &G_L) ;	  
      NbrFG = List_Nbr(G_L) ;
      FG = (int*)( G_L->array ) ;
      List_Read(FMMmat_S.Nd_L, i, &Nd_L) ;	  
      Nd = (int*)( Nd_L->array ) ;
      fprintf(file, "\nFAR (on %d) ", FMMmat_S.Obs);
      for (j = 0 ; j < NbrFG ; j++ )
       fprintf(file,"%d(%d) ", FG[j], Nd[j]);	
      fprintf(file, "\n\n");    
    }

  } /* iFMMEqu */
  
  fclose(file);
  
  GetDP_End;
}




void Print_FMM_Matrices( int DA, char* FileName ){
  
  int i, iDof, iEqu, iCom, iDir, iFMMEqu, NbrFMMEqu, FMMi ;
  int NbrGroups, NbrDir, NbrHar, NbrCom, NbrDof, NumDof, NbrEqu, NumEqu ;
  double **M ;

  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  List_T *NumDof_L, *NumEqu_L ;
  FILE * file ;
  

  GetDP_Begin("Print_FMM_Matrices");

  Msg(INFO, "Writing FMM '%s' ", FileName );
  
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  NbrDir = Current.FMM.NbrDir ; /* Current.FMM.N for Laplace3D */
  NbrHar = Current.NbrHar ;

  file = fopen(FileName, "w");  
  fprintf(file, "$FMM_Matrices \n\n");

  if (DA==0) fprintf(file, "$AGGREGATION \n"); else fprintf(file, "$DISAGGREGATION \n");
  fprintf(file, "NbrDirections = %d \n", NbrDir);
  fprintf(file, "NbrHarmonics = %d \n\n", NbrHar);

  for( iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;
    NbrCom = FMMmat_P->NbrCom ;
    if (DA == 0) FMMi = FMMmat_P->Src ; 
    else FMMi = FMMmat_P->Obs ;

    if (DA == 0) NbrGroups = List_Nbr( FMMmat_P->A_L ) ;
    else NbrGroups = List_Nbr( FMMmat_P->D_L ) ;

    fprintf(file, "NbrGroups (Support %d NbrComponents %d) = %d \n",
	    FMMi, NbrCom, NbrGroups);

    for (i = 0 ; i < NbrGroups ; i++){      
      if (DA==0){
	List_Read(FMMmat_P->A_L, i, &M ) ;
	List_Read(FMMmat_P->NumDof, i, &NumDof_L ) ;
	NbrDof = List_Nbr( NumDof_L) ;
	fprintf(file, "$matrixGROUP %d NbrDof %d\n", i, NbrDof);

	for (iDof = 0 ; iDof < NbrDof ; iDof++){
	  List_Read(NumDof_L, iDof, &NumDof);
	  fprintf(file, "$NumDof %d\n", NumDof);
	  
	  for (iDir = 0 ; iDir < NbrDir ; iDir++){
	    for (iCom = 0 ; iCom < NbrCom ; iCom++){ 
	      fprintf(file, "%d %3.8g %3.8g \n", iDir, 
		      M[iDof][iDir*NbrHar*NbrCom+iCom], M[iDof][iDir*NbrHar*NbrCom+NbrCom+iCom]);
	    }
	  }
	}   
      } 
      else {
	List_Read(FMMmat_P->D_L, i, &M ) ;
	List_Read(FMMmat_P->NumEqu, i, &NumEqu_L ) ;
	NbrEqu = List_Nbr(NumEqu_L) ;
	fprintf(file, "$matrixGROUP %d NbrEqu %d\n", i, NbrEqu);

	for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	  List_Read( NumEqu_L, iEqu, &NumEqu);
	  fprintf(file, "$NumEqu %d\n", NumEqu);
	  
	  for (iDir = 0 ; iDir < NbrDir ; iDir++){
	    for (iCom = 0 ; iCom < NbrCom ; iCom++){ 
	      fprintf(file, "%d %3.8g %3.8g \n", iDir, 
		      M[iEqu][iDir*NbrHar*NbrCom+iCom], M[iEqu][iDir*NbrHar*NbrCom+NbrCom+iCom]);
	    }
	  }
	}   
	
	
      }
      fprintf(file, "$END_matrixGROUP %d \n\n", i);
    } 
    fclose(file);
  }
  GetDP_End;
}



/* ------------------------------------------------------------------------- */ 
/* Print TRANSLATION matrix                                                  */ 
/* ------------------------------------------------------------------------- */ 

void Print_FMM_Translation( char* FileName ){

  int i, iFar, *FG, iDir, iDir2, iEqu ;
  int NbrEqu, NbrGroupSrc, NbrFG, NbrDir, NbrDir2, NbrHar ;
  double ***T;
  struct FMMGroup FMMGroup_S ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  List_T *FG_L ; 
  FILE * file ;


  GetDP_Begin("Print_FMM_Translation");

  Msg(INFO, "Writing FMM Translation matrix '%s' ", FileName );

  NbrDir = NbrDir2 = Current.FMM.NbrDir ;
  NbrHar = 2 ;
  NbrEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0) ;

  file = fopen(FileName, "w") ;
  fprintf(file, "$FMM_Translation Matrices \n\n") ;
  fprintf(file, "NbrDirections = %d \n", NbrDir) ;
  fprintf(file, "NbrHarmonics = %d \n\n", NbrHar) ;
 
  for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
    FMMmat_P = FMMmat_P0 + iEqu ;
    T = FMMmat_P->T ;

    fprintf(file, "SOURCE %d\n", FMMmat_P->Src) ; 
    fprintf(file, "OBSERV %d\n", FMMmat_P->Obs) ;

    List_Read(Problem_S.FMMGroup, FMMmat_P->Src, &FMMGroup_S);
    NbrGroupSrc = List_Nbr(FMMGroup_S.List) ;    
    fprintf(file, "NbrGroupSrc = %d \n", NbrGroupSrc);

    for (i = 0 ; i < NbrGroupSrc ; i++){
      List_Read(FMMmat_P->FG_L, i, &FG_L) ;	  
      NbrFG = List_Nbr(FG_L) ;
      FG = (int*)( FG_L->array ) ;
      
      for (iFar = 0 ; iFar < NbrFG ; iFar++){
	fprintf(file, "$GROUPS Origin %d (Support %d) Destination %d (Support %d) \n",
		i, FMMmat_P->Src, FG[iFar], FMMmat_P->Obs);      
	for (iDir = 0 ; iDir < NbrDir ; iDir++)
	  for (iDir2 = 0 ; iDir2 < NbrDir2 ; iDir2++)	
	    fprintf(file, "%d %3.8g %3.8g \n",
		    iDir, T[i][iFar][iDir*NbrHar*NbrDir+iDir2*NbrHar ], T[i][iFar][iDir*NbrHar*NbrDir+iDir2*NbrHar+1]) ;
      }      
      fprintf(file, "$END_GROUP Origin %d \n\n", i);
    }
  }

  fclose(file);
    
  GetDP_End;
}



/* ------------------------------------------------------------------------- */ 
/* Print FMM DTA matrix                                                      */ 
/* ------------------------------------------------------------------------- */ 

void Print_FMM_MatrixDTA( int n, double **DTA, char* FileName ){

  int  iDof, iEqu ; 
  FILE * file ;

  GetDP_Begin("Print_FMM_MatrixDTA");

  Msg(INFO, "Writing FMM DTA matrix : '%s' ", FileName );

  file = fopen(FileName, "w");
  
  fprintf(file, "$FMM_Matrix_DTA  NbrDof = %d\n", n) ;
  
  /* GetDP works with the transposed system matrix */
  for (iEqu = 0 ; iEqu < n ; iEqu++) 
    for (iDof = 0 ; iDof < n ; iDof++)   
      fprintf(file, "%d %d %.8g\n", iEqu+1, iDof+1, DTA[iEqu][iDof]) ; 
  
  fclose(file);

  GetDP_End;
}