print_fmm.c

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

#define RCSID "$Id: Print_FMM.c,v 1.12 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 "GeoData.h"
#include "Numeric.h"
#include "Data_Passive.h"
#include "CurrentData.h"
#include "Tools.h"
#include "ExtendedGroup.h"
#include "Data_FMM.h"

int Get_GmshElementType(int Type) ;

void  Geo_WriteFileFMMGroups(struct GeoData *GeoData_P, char * FileName){

  int  i, j, i_Support, NbrInSupport, NbrGroups, NbrElmsGroupi, k = -1 ;
  int *ElmsGroupi_P0 ;
  struct Geo_Element *Element_P ;
  struct FMMData  *FMMData_P0 ;
  struct FMMGroup  FMMGroup_S ;

  FILE * file ;

  GetDP_Begin("Geo_WriteFileFMMGroups");

  Msg(INFO, "Writing FMMGroups File .pos : '%s' ", FileName );
  file = fopen(FileName, "w");

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

    fprintf(file, "View \"FMM Support %d (%d)\" { \n", i_Support, NbrGroups ) ;
    
    for( i = 0 ; i < NbrGroups ; i++ ){
      if ( k < 10 ) k++;
      else k = 0 ;
      NbrElmsGroupi =  List_Nbr((FMMData_P0+i)->Element) ;
      ElmsGroupi_P0 = (int*)List_Pointer((FMMData_P0+i)->Element,0) ;
      for( j = 0 ; j < NbrElmsGroupi ; j++){
	Element_P = Geo_GetGeoElementOfNum(ElmsGroupi_P0[j]);
	switch(Element_P->Type) {
	case POINT : 
	  fprintf(file, "SP( %f , %f , %f ) { %d };\n",
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->z, i); break;	  
	case LINE :
	  fprintf(file, "SL( %f , %f , %f , %f , %f , %f ) { %d ,%d };\n",
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->z,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->z, i, i); break;
	case TRIANGLE : 
	  fprintf(file, "ST( %f , %f , %f , %f , %f , %f, %f , %f , %f ) { %d, %d ,%d };\n",
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->z,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->z,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[2])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[2])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[2])->z, k, k, k); break;
	case QUADRANGLE : 
	  fprintf(file, "SQ( %f , %f , %f , %f , %f , %f, %f , %f , %f, %f , %f , %f) { %d, %d ,%d, %d };\n",
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[0])->z,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[1])->z,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[2])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[2])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[2])->z,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[3])->x,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[3])->y,
		  Geo_GetGeoNodeOfNum(Element_P->NumNodes[3])->z, i, i, i, i); break;
	default :  Msg(GERROR, "Wrong Type of Element for FMM groups"); break;   
	}
      }
    } 
    fprintf(file, "};\n");    
  }

  fclose(file);

  GetDP_End ;
}
  


void  Geo_WriteFileFMMGroupsCenter( char * FileName ){

  int NbrInSupport, i_Support, NbrGroups, i ;
  struct FMMData *FMMData_P0 ;
  struct FMMGroup  FMMGroup_S ;
  List_T *FMMGroup_L ;  
 
  FILE * file ;


  GetDP_Begin("Geo_WriteFileFMMGroupsCenter");

  Msg(INFO, "Writing FMMGroupsCenters File .pos : '%s' ", FileName );
  file = fopen(FileName, "w");

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

    fprintf(file, "View \"FMMCenters Support %d \" { \n", i_Support);

    for (i=0 ; i < NbrGroups ; i++)
      fprintf(file, "SP( %f , %f , %f ) { %d };\n",
	      (FMMData_P0+i)->Xgc, (FMMData_P0+i)->Ygc, (FMMData_P0+i)->Zgc, i);
  
    fprintf(file, "};\n");
  }    
  fclose(file);
  
  GetDP_End ;
}



void  Geo_WriteFileMshFMMGroups( struct GeoData *GeoData_P, char * FileName ){

  int i, j, iNd, i_Support, NbrElmsInSupport, NbrNdsInSupport ;
  int NbrInSupport, NbrNds, NbrGroups, NbrElmsGroupi ;
  int ElemType, *ElmsGroupi_P0 ;
  struct Geo_Node *Nodes_P0 ;
  struct Geo_Element *Element_P ;
  struct FMMData *FMMData_P0 ;
  struct FMMGroup  FMMGroup_S ;
  List_T *NodesInRegion ;
  List_T *FMMGroup_L ;  
  FILE * file ;

  GetDP_Begin("Geo_WriteFileMshFMMGroups");

  Msg(INFO, "Writing FMMGroups  File .msh : '%s' ", FileName );
  NbrNds = List_Nbr(GeoData_P->Nodes) ;
  Nodes_P0 = (struct Geo_Node*)List_Pointer( GeoData_P->Nodes, 0 ) ;  


  file = fopen(FileName, "w") ;

  NbrInSupport = List_Nbr(Problem_S.FMMGroup) ;

  for(i_Support = 0 ; i_Support < NbrInSupport ; i_Support++ ){
    NbrElmsInSupport = 0 ; 
    List_Read(Problem_S.FMMGroup, i_Support, &FMMGroup_S) ;
    FMMGroup_L = FMMGroup_S.List ;
   
    Generate_ElementaryEntities( ((struct Group *) List_Pointer(Problem_S.Group, FMMGroup_S.InIndex))->InitialList ,
				 &NodesInRegion, NODESOF) ;
    NbrNdsInSupport = List_Nbr(NodesInRegion);
    
    fprintf(file, "$NOD \n");
    fprintf(file, "%d \n", NbrNdsInSupport);    
    for (iNd = 0 ; iNd < NbrNds ; iNd++)
      if(List_Search(NodesInRegion, &(Nodes_P0+iNd)->Num, fcmp_int))
	fprintf(file, "%d %f %f %f \n", (Nodes_P0+iNd)->Num,(Nodes_P0+iNd)->x, (Nodes_P0+iNd)->y, (Nodes_P0+iNd)->z);
    fprintf(file, "$ENDNOD \n");

    NbrGroups = List_Nbr(FMMGroup_L) ;
    FMMData_P0 = (struct FMMData*)List_Pointer(FMMGroup_L,0) ;  

    for( i = 0 ; i < NbrGroups ; i++)
      NbrElmsInSupport += List_Nbr((FMMData_P0+i)->Element) ;

    fprintf(file, "$ELM Support %d \n", i_Support);
    fprintf(file, "%d \n", NbrElmsInSupport);
 
    for( i = 0 ; i < NbrGroups ; i++){
      NbrElmsGroupi =  List_Nbr((FMMData_P0+i)->Element) ;
      ElmsGroupi_P0 = (int*)List_Pointer((FMMData_P0+i)->Element,0) ;
      
      for (j=0 ; j < NbrElmsGroupi ; j++){
	Element_P = Geo_GetGeoElementOfNum(ElmsGroupi_P0[j]);
	ElemType = Get_GmshElementType(Element_P->Type);
	switch(Element_P->Type){
	case POINT :
	  fprintf(file, "%d %d %d %d %d %d \n", 
		  Element_P->Num , ElemType, i, Element_P->Region, Element_P->NbrNodes,
		  Element_P->NumNodes[0]); break;
	case LINE :
	  fprintf(file, "%d %d %d %d %d %d %d \n", 
		  Element_P->Num , ElemType, i, Element_P->Region, Element_P->NbrNodes,
		  Element_P->NumNodes[0], Element_P->NumNodes[1]); break;
	case TRIANGLE : 
	  fprintf(file, "%d %d %d %d %d %d %d %d \n", 
		  Element_P->Num , ElemType, i, Element_P->Region, Element_P->NbrNodes,
		  Element_P->NumNodes[0], Element_P->NumNodes[1] ,Element_P->NumNodes[2]); break;
	case QUADRANGLE : 
	  fprintf(file, "%d %d %d %d %d %d %d %d %d \n", 
		  Element_P->Num , ElemType, i, Element_P->Region, Element_P->NbrNodes,
		  Element_P->NumNodes[0], Element_P->NumNodes[1],
		  Element_P->NumNodes[2],Element_P->NumNodes[3]); break;

	default :  Msg(GERROR, "Wrong Type of Element for FMM groups"); break;
	}
      }
    }
    fprintf(file, "$ENDELM\n");
  }
  
  fclose(file);

  GetDP_End;
}



void Print_FMMGroupInfo( char* FileName ){
  
  int i, j, NbrInSupport, NumElm, *NG, *FG, *Nd, NbrGroupsSrc, NbrGroupsObs ;
  int NbrDir, NbrHar, NbrCom, NbrElmsGroup, NbrNG, NbrFG ;
  int NbrDof=0, NbrEqu=0, *NumDof, *NumEqu, NbrFMMEqu, iFMMEqu ;
  double  StorageFMM=0., StorageFMMA, StorageFMMD, StorageFMMT, StorageDir=0., TotalMem ;
  int NbrMultLaplace2D, NbrMultLaplace3D, N ;
  int NbrDirMAX = 0, NbrDirAV = 0 ;
  struct FMMData *FMMDataObs_P0, *FMMDataSrc_P0 ;
  struct FMMGroup FMMGroupObs_S, FMMGroupSrc_S ;
  struct FMMmat FMMmat_S ;
  List_T *G_L, *Nd_L, *NumDof_L, *NumEqu_L ;
 

  FILE * file ;
  

  GetDP_Begin("Print_FMMGroupInfo") ;
  
  Msg(INFO, "Writing FMM information : '%s' ", FileName );
 
  NbrDir = Current.FMM.NbrDir ;
  N = (NbrDir-1)*(NbrDir+1);
  NbrHar = 2 ;

  if (Current.FMM.Precision == 0) NbrDirMAX = NbrDirAV = NbrDir ;

  NbrCom = Current.FMM.NbrCom ;
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  NbrInSupport = List_Nbr(Problem_S.FMMGroup) ;

  NbrMultLaplace2D = NbrMultLaplace3D = 0 ;

  file = fopen(FileName, "w");
  fprintf(file, "$FMMGroupsInfo\n");
  fprintf(file, "NbrDirections = %d NbrHarmonics = %d NbrCom = %d (Warning: Default == SCALAR)\n", NbrDir, NbrHar, NbrCom) ;
  fprintf(file, "NbrFMMEqu = %d NbrSupports = %d\n\n", NbrFMMEqu, NbrInSupport) ;

  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) ;
  }

  TotalMem = 0. ; 

  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu++){
    StorageFMM = StorageDir = 0. ;
    StorageFMMA = StorageFMMD = StorageFMMT = 0. ;

    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) ;
      fprintf(file, "FmmGROUP #Elm #Dof #Equ\n");
     
      for (i = 0 ; i < NbrGroupsSrc ; i++){
	NbrElmsGroup =  List_Nbr((FMMDataSrc_P0+i)->Element) ;
	List_Read(FMMmat_S.NumDof, i, &NumDof_L) ;
	NbrDof = List_Nbr(NumDof_L) ;
	List_Read(FMMmat_S.NumEqu, i, &NumEqu_L) ;
	NbrEqu = List_Nbr(NumEqu_L) ;
	fprintf(file, "%5d  %5d  %5d  %5d\n", i, NbrElmsGroup, NbrDof, NbrEqu);
	StorageFMMA += NbrDof*NbrCom*NbrHar*NbrDir* 32./1024. ;
	StorageFMMD += NbrEqu*NbrCom*NbrHar*NbrDir* 32./1024. ;

	NbrMultLaplace2D += NbrDof*NbrDir*NbrHar + NbrEqu *NbrDir *NbrHar + NbrDir*NbrDir*6 ;
	NbrMultLaplace3D += NbrDof*N*NbrHar + NbrEqu * N *NbrHar + N*N*6 ;
	
	List_Read(FMMmat_S.FG_L, i, &G_L) ;
       
	NbrFG = List_Nbr(G_L) ;
	StorageFMMT += NbrFG*NbrHar * N * 32./1024. ;
      }
    }   
    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, "FmmGROUP   #Elm  #Dof\n");
      for (i = 0 ; i < NbrGroupsSrc ; i++){
	NbrElmsGroup =  List_Nbr((FMMDataSrc_P0+i)->Element) ;
	List_Read(FMMmat_S.NumDof, i, &NumDof_L) ;
	NbrDof = List_Nbr(NumDof_L) ;
	fprintf(file, "%5d  %5d  %5d\n", i, NbrElmsGroup, NbrDof);
	StorageFMMA += NbrDof*NbrCom*NbrHar*NbrDir* 32./1024. ;
	List_Read(FMMmat_S.FG_L, i, &G_L) ;	  
	NbrFG = List_Nbr(G_L) ;
	StorageFMMT += NbrFG*NbrHar*NbrDir* 32./1024. ;

      }
      fprintf(file, "FmmGROUP    #Elm  #Equ\n");
      for (i = 0 ; i < NbrGroupsObs ; i++){
	NbrElmsGroup =  List_Nbr((FMMDataObs_P0+i)->Element) ;
	List_Read(FMMmat_S.NumEqu, i, &NumEqu_L) ;
	NbrDof = List_Nbr(NumEqu_L) ;
	fprintf(file, "%5d  %5d  %5d\n", i, NbrElmsGroup, NbrEqu);
	StorageFMMD += NbrEqu*NbrCom*NbrHar*NbrDir* 32./1024. ;
      }
    }

    StorageFMM += StorageFMMA + StorageFMMD + StorageFMMT ;   
  
  

    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) ;
      StorageDir += (NbrDof+NbrEqu)*(NbrDof+NbrEqu)*NbrHar*32./1024. ; 
      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 ) ;
      fprintf(file, "\nFAR (on %d) ", FMMmat_S.Obs);
     
      if(Current.FMM.Precision != 0){
	List_Read(FMMmat_S.Nd_L, i, &Nd_L) ;	  
	Nd = (int*)( Nd_L->array ) ;      
	for (j = 0 ; j < NbrFG ; j++ ){
	  fprintf(file,"%d(%d) ", FG[j], Nd[j]);	
	  NbrDirAV += Nd[j] ;
	  NbrDirMAX = (NbrDirMAX < Nd[j]) ? Nd[j] : NbrDirMAX ; 
	}
      }
      else 
	for (j = 0 ; j < NbrFG ; j++ )