#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 <getdp@geuz.org>.
 *
 * 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++ )
	  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;
}