#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 <getdp@geuz.org>.
 *
 * 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);
      }
      List_Reset(FMMmat_P->A_L);
      List_Reset(FMMmat_P->D_L);
    }

  }
  
  for (i_FMMEqu = 0 ; i_FMMEqu < NbrFMMEqu ; i_FMMEqu++){
    FMMmat_P  = FMMmat_P0 + i_FMMEqu ;  
    i_Obs = FMMmat_P->Obs ;
    i_Src = FMMmat_P->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) ;
      }   
      
      if(Flag_FMMDA){
	for(i = 0 ; i < NbrGroupSrc ; i++) List_Add(FMMmat_P->A_L, &Mat) ;
	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) ;
      }

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

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

	List_Add( FMMmat_P->NG_L, &NGi );
	List_Add( FMMmat_P->FG_L, &FGi );
	List_Add( FMMmat_P->Nd_L, &Ndi );
      }/* NbrGroupSrc */
    }
  }/* i_FMMEqu */

  GetDP_End;

}


void  Geo_SetXYZMinMaxInElement(int Num, double *Xmin, double *Ymin, double *Zmin,
				double *Xmax, double *Ymax, double *Zmax){
 
  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_SetXYZMinMaxInElement");

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

  for (j = 0;  j < Element->NbrNodes ; j++){
    *Xmax = (*Xmax >= x[j]) ? *Xmax : x[j] ;
    *Ymax = (*Ymax >= y[j]) ? *Ymax : y[j] ;
    *Zmax = (*Zmax >= z[j]) ? *Zmax : z[j] ;
    *Xmin = (*Xmin <= x[j]) ? *Xmin : x[j] ;
    *Ymin = (*Ymin <= y[j]) ? *Ymin : y[j] ;
    *Zmin = (*Zmin <= z[j]) ? *Zmin : z[j] ;
  }
  
  GetDP_End ;
}



double Geo_SetMaxEdgeLength(int Num){
 
  int j ;
  double D = 0. ,Dmax =0. ;
  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_SetMaxEdgeLength");

  Element =  Geo_GetGeoElementOfNum(Num) ;    
  Geo_GetNodesCoordinates(Element->NbrNodes, Element->NumNodes,
			  x, y, z);
  
  for (j = 0;  j < Element->NbrNodes ; j++){
    if (j+1<Element->NbrNodes)
      D = sqrt(SQU(x[j]-x[j+1])+ SQU(y[j]-y[j+1])+ SQU( z[j]-z[j+1])) ;
  
    Dmax = ( D > Dmax ) ? D : Dmax ;
  }
  
  GetDP_Return(Dmax) ;
}


double Geo_GetMaxDistance2Elms(int Numi, int Numj){
 
  int i, j ;
  double D = 0. ,Dmax =0. ;
  double x1[NBR_MAX_NODES_IN_ELEMENT], y1[NBR_MAX_NODES_IN_ELEMENT], z1[NBR_MAX_NODES_IN_ELEMENT];
  double x2[NBR_MAX_NODES_IN_ELEMENT], y2[NBR_MAX_NODES_IN_ELEMENT], z2[NBR_MAX_NODES_IN_ELEMENT];

  struct Geo_Element *Elementi, *Elementj ;  

  GetDP_Begin("Geo_GetMaxDistance2Elms");

  Elementi =  Geo_GetGeoElementOfNum(Numi) ;    
  Geo_GetNodesCoordinates(Elementi->NbrNodes, Elementi->NumNodes,
			  x1, y1, z1);

  Elementj =  Geo_GetGeoElementOfNum(Numj) ;    
  Geo_GetNodesCoordinates(Elementj->NbrNodes, Elementj->NumNodes,
			  x2, y2, z2);
  
  for (i = 0;  i < Elementi->NbrNodes ; i++)
    for (j = 0;  j < Elementj->NbrNodes ; j++){
      D = sqrt(SQU(x1[i]-x2[j])+ SQU(y1[i]-y2[j])+ SQU(z1[i]-z2[j])) ;  
      Dmax = ( D > Dmax ) ? D : Dmax ;
    }
  
  GetDP_Return(Dmax) ;
}



int FMM_NbrSpectralDirections( double k0 ){

  int i, j, i_Support, NbrInSupport, NbrGroups, NbrElmsGroupi ;
  int NumElmj, NumElmjp ;
  double D = 0., Daux = 0., Centroidj[3], Centroidjp[3] ; 
  struct FMMData  *FMMData_P0 ;
  struct FMMGroup  FMMGroup_S ;

  GetDP_Begin("FMM_NbrSpectralDirections");

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

    for(i = 0; i < NbrGroups; i++ ){
      NbrElmsGroupi = List_Nbr((FMMData_P0+i)->Element) ;
      
      if (NbrElmsGroupi==1){ 
	List_Read((FMMData_P0+i)->Element, 0, &NumElmj);
	Daux = Geo_SetMaxEdgeLength(NumElmj);
      }
      else{
	for (j = 0 ; j < NbrElmsGroupi ; j++){
	  if ((j+1) < NbrElmsGroupi){
	    List_Read((FMMData_P0+i)->Element, j, &NumElmj);
	    List_Read((FMMData_P0+i)->Element, j+1, &NumElmjp);
	    
	    Geo_SetCentroidCoordinates(NumElmj, Centroidj);
	    Geo_SetCentroidCoordinates(NumElmjp, Centroidjp);
	    
	    Daux = sqrt(SQU(Centroidj[0]-Centroidjp[0])+ SQU(Centroidj[1]-Centroidjp[1])+ SQU(Centroidj[2]-Centroidjp[2])) ;       
	  }
	  D = (Daux>D) ? Daux : D ;
	}
      }
      D = (Daux>D) ? Daux : D ;
    }    
  }
  
  Msg(INFO,"D_FMM = %.8f", D);
  
  GetDP_Return( (int)ceil( k0*D+5.*log(k0*D+PI) ) );
}




int FMM_SetTruncation( double Rsrc, double Robs, double D, int Dimension){

  int NbrDir=0, i, j, k, exp ;
  double Precision ; 

  GetDP_Begin("FMM_SetTruncation");

  Precision = Current.FMM.Precision ;
  exp = - (int)floor(log10(Precision)) ;

  i = (int)floor(Rsrc/D/0.05)-1 ; 
  j = (int)floor(Robs/D/0.05)-1 ;
  if (i<0) i = 0 ;
  if (j<0) j = 0 ;
  k = ((i + j * 7)>48)?48: i+j*7 ;


  if (Dimension == _2D){
    switch(exp){
    case 3 :  NbrDir = t2D1e_3[k] ; break ;  
    case 4 :  NbrDir = t2D1e_4[k] ; break ;  
    case 5 :  NbrDir = t2D1e_5[k] ; break ;  
    case 6 :  NbrDir = t2D1e_6[k] ; break ;  
    case 7 :  NbrDir = t2D1e_7[k] ; break ;  
    case 8 :  NbrDir = t2D1e_8[k] ; break ;  
    case 9 :  NbrDir = t2D1e_9[k] ; break ;
    default: 
      Msg(GERROR, "Precision for FMM truncation can only be: 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9");
    }
    if (NbrDir > 11){ if (Flag_FMM == 1) Msg(WARNING,"NbrDir = %d changed to 11 for memory requirements", NbrDir) ; NbrDir = 11 ;}
  }
  else{
    switch(exp){
    case 3 :  NbrDir = t3D1e_3[k] ; break ;  
    case 4 :  NbrDir = t3D1e_4[k] ; break ;  
    case 5 :  NbrDir = t3D1e_5[k] ; break ;  
    case 6 :  NbrDir = t3D1e_6[k] ; break ;  
    case 7 :  NbrDir = t3D1e_7[k] ; break ;  
    case 8 :  NbrDir = t3D1e_8[k] ; break ;  
    case 9 :  NbrDir = t3D1e_9[k] ; break ;
    default: 
      Msg(GERROR, "Precision for FMM truncation can only be: 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9");
  }
    if (NbrDir > 15){ 
      if (Flag_FMM == 1) 
	Msg(WARNING,"NbrDir = %d changed to 15 for memory requirements Rsrc/D %g Robs/D %g", NbrDir,Rsrc/D,Robs/D) ; NbrDir = 15 ;}
  }
 

  if (Flag_FMM == 2)  
    NbrDir = (NbrDir > Current.FMM.NbrDir) ? Current.FMM.NbrDir : NbrDir ; 
  else
    Current.FMM.NbrDir = (NbrDir>Current.FMM.NbrDir) ? NbrDir : Current.FMM.NbrDir ;
  
  GetDP_Return(NbrDir);
}


int FMM_SetTruncationOLD( double Rmax, double D, int Dimension){

  int NbrDir ;
  double Precision ; 

  GetDP_Begin("FMM_SetTruncation");

 /* Rmax == Radius of the circle/sphere enclosing the FMM group */  
  Precision = Current.FMM.Precision ;

  if (Dimension == _2D){
    NbrDir =  (int)floor(-log(Precision)/log(D/Rmax)) ; /* Truncation parameter */
    if (NbrDir > 15){ if (Flag_FMM == 1) Msg(WARNING,"NbrDir = %d changed to 15 for memory requirements", NbrDir) ; NbrDir = 15 ;}
  }
  else{
    NbrDir =  (int)floor(-log(Precision)/log(2*D/Rmax)) ;
    if (NbrDir > 20){ if (Flag_FMM == 1) Msg(WARNING,"NbrDir = %d changed to 20 for memory requirements", NbrDir) ; NbrDir = 20 ;}
  }
 

  if (Flag_FMM == 2)  
    NbrDir = (NbrDir > Current.FMM.NbrDir) ? Current.FMM.NbrDir : NbrDir ; 
  else
    Current.FMM.NbrDir = (NbrDir>Current.FMM.NbrDir) ? NbrDir : Current.FMM.NbrDir ;
  
  GetDP_Return(NbrDir);
}



int  NeighbouringGroups( int FMMGroupEObs, int FMMGroupESrc ){

  int iFMMEqu, NbrFMMEqu, FMMObs, FMMSrc, NbrNG=0, i, *NG ;
  struct FMMmat  *FMMmat_P0 ;
  List_T *NG_L ;

  GetDP_Begin("NeigbouringGroups");

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0) ;
  
  iFMMEqu = 0 ; 
  while (iFMMEqu<NbrFMMEqu && ((FMMmat_P0+iFMMEqu)->Obs != Current.FMM.Obs || (FMMmat_P0+iFMMEqu)->Src != Current.FMM.Src)) 
    iFMMEqu++ ;

  if ( iFMMEqu == NbrFMMEqu )  Msg(GERROR, "Wrong Supports for FMM: Source %d Observation %d", Current.FMM.Src, Current.FMM.Obs) ; 

  FMMObs = (FMMmat_P0+iFMMEqu)->Obs ;
  FMMSrc = (FMMmat_P0+iFMMEqu)->Src ;

  if (FMMObs == FMMSrc &&  FMMGroupEObs == FMMGroupESrc) GetDP_Return(1) ;


  List_Read((FMMmat_P0+iFMMEqu)->NG_L, FMMGroupESrc, &NG_L);
  NbrNG = List_Nbr(NG_L) ;  

  NG = (int*)(NG_L->array) ;   
  
  for ( i = 0 ; i < NbrNG ; i++ )
    if (NG[i] == FMMGroupEObs) GetDP_Return(1) ; 
     
  GetDP_Return(0) ;
}



int  Get_NextElementSourceNeighbour( struct Element      * Element ) {

  struct Element * ElementSource ;

  GetDP_Begin("Get_NextElementSourceNeighbour");
  
  ElementSource = Element->ElementSource ;

  while (Get_NextElementSource(ElementSource))
    if (NeighbouringGroups(Element->FMMGroup, ElementSource->FMMGroup)) GetDP_Return(1) ;
  
  GetDP_Return(0) ;
}


int  Get_NextElementSourceInGroup(struct Element *Element) {
  struct Element * ElementSource ;

  GetDP_Begin("Get_NextElementSourceInGroup");

  ElementSource = Element->ElementSource ;
  while (Get_NextElementSource(ElementSource))
    if (ElementSource->FMMGroup == Element->FMMGroup) GetDP_Return(1) ;  
  GetDP_Return(0) ;
}



void  Geo_CreateFMMGroup( int InSupport, struct GeoData *GeoData_P, double k0 ){

  int NbrGroupsInSupport, iGroup, NbrFMMGroupsInRegion ;
  int NbrElms, NbrElmsInRegion, NbrElmsInSupport ;
  int j, k, l, m, FoundHome, iFMMGroup, iElm ;
  int ***NbrElmsDiv, ****ElmsDiv, Dimension ;
  double  Xmax, Xmin, Ymax, Ymin, Zmax, Zmin, l0 ;
  int NbrDivX = 1, NbrDivY = 1, NbrDivZ = 1 ;
  double  x1, x2, y1, y2, z1, z2, Centroid[3], tol ;
  struct Value DivXYZ ;
  struct Geo_Element *Elements_P0 ;
  struct FMMData FMMData_S ;
  struct FMMGroup FMMGroup_S ;
  List_T * FMMGroup_L ;  


  GetDP_Begin("Geo_CreateFMMGroup");

  /* Dimension = (int)GeoData_P->Dimension ; */
  Dimension = Current.FMM.Dimension ;    

  tol = 1e-6;

  if (!List_Nbr(Problem_S.FMMGroup)){
    if(k0 > 0){
      l0 = k0/TWO_PI ; /* Wavelength */
      Msg(INFO,"Helmholtz %dD, Wavelength = %g", Dimension, l0);
    }  
    else{
      Msg(INFO,"Laplace %dD", Dimension);
    }
  }

  NbrElms = List_Nbr(GeoData_P->Elements) ;
  Elements_P0 =  (struct Geo_Element*)List_Pointer(GeoData_P->Elements, 0) ;

  FMMGroup_L = List_Create (10, 10, sizeof(struct FMMData)) ; 
  
  NbrElmsInSupport = 0 ;  
  NbrGroupsInSupport = List_Nbr(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList ); 

  iFMMGroup = 0 ;
  
  for (iGroup = 0 ; iGroup < NbrGroupsInSupport ; iGroup++){
    List_Read(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList, iGroup, &Current.Region) ;    
    Get_ValueOfExpressionByIndex(Current.FMM.DivXYZIndex, NULL, 0., 0., 0., &DivXYZ) ;
    NbrDivX = (int)DivXYZ.Val[0] ;
    NbrDivY = (int)DivXYZ.Val[1] ;
    NbrDivZ = (int)DivXYZ.Val[2] ;
  
    Msg(INFO,"Region = %d DivX = %d  DivY = %d  DivZ = %d", Current.Region, NbrDivX, NbrDivY, NbrDivZ) ; 
    
    NbrElmsInRegion = 0 ; 
    Xmax = Ymax = Zmax = Xmin = Ymin = Zmin = 0. ;

    for(iElm = 0 ; iElm < NbrElms ; iElm++){
      if((Elements_P0+iElm)->Region == Current.Region){
	NbrElmsInRegion++ ; 
	Geo_SetXYZMinMaxInElement((Elements_P0+iElm)->Num, &Xmin, &Ymin, &Zmin, &Xmax, &Ymax, &Zmax ) ;      
      }      
    }

    NbrElmsInSupport += NbrElmsInRegion ; 

    NbrElmsDiv = (int ***)Malloc(NbrDivX*sizeof(int**));
    for (j = 0 ; j < NbrDivX ; j++){
      NbrElmsDiv[j] = (int **)Calloc(NbrDivY, sizeof(int*));
      for (k = 0 ; k < NbrDivY ; k++)
	NbrElmsDiv[j][k] = (int *)Calloc(NbrDivZ, sizeof(int));
    }
    
    ElmsDiv = (int ****)Malloc(NbrDivX*sizeof(int***));
    for (j = 0 ; j < NbrDivX ; j++){
      ElmsDiv[j] = (int ***)Malloc(NbrDivY*sizeof(int**)); 
      for (k = 0 ; k < NbrDivY ; k++){
	ElmsDiv[j][k] = (int **)Malloc(NbrDivZ*sizeof(int*));
	for (l = 0 ; l < NbrDivZ ; l++)
	  ElmsDiv[j][k][l] = (int *)Malloc(sizeof(int));
      }
    }

    for (iElm = 0 ; iElm < NbrElms ; iElm++){ 
      FoundHome = 0;
      
      for (j = 0 ; j < NbrDivX ; j++){
	for (k = 0 ; k < NbrDivY ; k++){
	  for (l = 0 ; l < NbrDivZ ; l++){
	    x1 = -tol+Xmin +(double)j*(Xmax-Xmin)/(double)NbrDivX;
	    x2 =  tol+Xmin +(double)(j+1)*(Xmax-Xmin)/(double)NbrDivX;
	    y1 = -tol+Ymin +(double)k*(Ymax-Ymin)/(double)NbrDivY;
	    y2 =  tol+Ymin +(double)(k+1)*(Ymax-Ymin)/(double)NbrDivY;
	    z1 = -tol+Zmin +(double)l*(Zmax-Zmin)/(double)NbrDivZ;
	    z2 =  tol+Zmin +(double)(l+1)*(Zmax-Zmin)/(double)NbrDivZ;

	    if( Current.Region == (Elements_P0+iElm)->Region ){
	      Geo_SetCentroidCoordinates((Elements_P0+iElm)->Num, Centroid);
	      
	      if (x1 <= Centroid[0] && Centroid[0] <= x2 &&
		  y1 <= Centroid[1] && Centroid[1] <= y2 &&
		  z1 <= Centroid[2] && Centroid[2] <= z2 ){
		if(!NbrElmsDiv[j][k][l]){
		  NbrElmsDiv[j][k][l] = 1;
		  ElmsDiv[j][k][l][0] = iElm ;
		} else {
		  NbrElmsDiv[j][k][l] += 1;
		  ElmsDiv[j][k][l] = (int *)realloc(ElmsDiv[j][k][l], NbrElmsDiv[j][k][l]*sizeof(int));
		  ElmsDiv[j][k][l][NbrElmsDiv[j][k][l]-1] = iElm;	    
		}
		FoundHome = 1; break;
	      }
	    }/* if Region */
	    if(FoundHome)break;
	  }
	  if(FoundHome)break;
	}
	if(FoundHome)break;
      }
      if(!FoundHome && (Current.Region == (Elements_P0+iElm)->Region)) 
	Msg(WARNING,"Element %d has not been assigned to a FMMGroup", iElm);
    }/* Elms NbrElms */ 

    NbrFMMGroupsInRegion = 0;
    for (j = 0 ; j < NbrDivX ; j++)
      for (k = 0 ; k < NbrDivY ; k++)
	for (l = 0 ; l < NbrDivZ ; l++)
	  if(NbrElmsDiv[j][k][l]){
	    Geo_InitFMMData(&FMMData_S) ;
	    FMMData_S.Element = List_Create(NbrElmsDiv[j][k][l], 1, sizeof(int)) ;
	    for(m = 0 ; m < NbrElmsDiv[j][k][l] ; m++){
	      List_Add(FMMData_S.Element, &(Elements_P0 + ElmsDiv[j][k][l][m])->Num) ;
	      Geo_SetCentroidCoordinates((Elements_P0+ElmsDiv[j][k][l][m])->Num, Centroid) ;

	      FMMData_S.Xgc += Centroid[0] ;
	      FMMData_S.Ygc += Centroid[1] ;
	      FMMData_S.Zgc += Centroid[2] ;
	      
	      (Elements_P0 + ElmsDiv[j][k][l][m])->FMMGroup = iFMMGroup ;
	    }	
	    
	    FMMData_S.Group = iFMMGroup;  
	    FMMData_S.Xgc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Ygc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Zgc /= NbrElmsDiv[j][k][l];
	    List_Add(FMMGroup_L, &FMMData_S) ;
	    
	    iFMMGroup +=1 ;
	    NbrFMMGroupsInRegion += 1 ;
	  }/* if( NbrElmsDiv[j][k][l] ) */    
  
    
    Free(NbrElmsDiv);
    Free(ElmsDiv);

    Msg(INFO,"SupportIndex = %d  Region = %d  NbrElmsInRegion = %d  NbrFMMGroupsInRegion = %d",
	InSupport, Current.Region, NbrElmsInRegion, NbrFMMGroupsInRegion) ; 
  }/* iGroup NbrGroupInSuppport*/
 
  Msg(INFO,"SupportIndex = %d  NbrElmsInSupport = %d  NbrFMMGroupsInSupport = %d",
      InSupport, NbrElmsInSupport, iFMMGroup) ; 

  FMMGroup_S.InIndex = InSupport ;
  FMMGroup_S.List = FMMGroup_L ;
  List_Add(Problem_S.FMMGroup, &FMMGroup_S ); 


  GetDP_End ;
}


#if 0 /* FIXME: this is just a sketch -- is not meant to work or even compile */
  
void  Geo_CreateMultilevelFMMGroup( int InSupport, struct GeoData *GeoData_P, double k0 ){

  int NbrGroupsInSupport, iGroup, NbrFMMGroupsInRegion ;
  int NbrElms, NbrElmsInRegion, NbrElmsInSupport ;
  int j, k, l, m, FoundHome, iFMMGroup, iElm, iFMMParentGroup ;
  int ***NbrElmsDiv, ****ElmsDiv, ***NbrGrsDiv, ****GrsDiv ;
  double  Xmax, Xmin, Ymax, Ymin, Zmax, Zmin ;
  int NbrDivX, NbrDivY, NbrDivZ;
  int NbrDivX0 = 1, NbrDivY0 = 1, NbrDivZ0 = 1, NbrDivX1 = 1, NbrDivY1 = 1, NbrDivZ1 = 1  ;
  double  x1, x2, y1, y2, z1, z2, Centroid[3], tol ;
  struct Value DivXYZ ;
  struct Geo_Element *Elements_P0 ;
  struct FMMData FMMData_S ;
  struct FMMData FMMParentData_S ;
  struct FMMGroup FMMGroup_S ;

  List_T * FMMGroup_L ;  
  List_T * FMMParentGroup_L ;  


  GetDP_Begin("Geo_CreateMultilevelFMMGroup");

  tol = 1e-6;

  NbrElms = List_Nbr(GeoData_P->Elements) ;
  Elements_P0 =  (struct Geo_Element*)List_Pointer(GeoData_P->Elements, 0) ;

  FMMGroup_L = List_Create (10, 10, sizeof(struct FMMData)) ; 
  FMMParentGroup_L = List_Create (10, 10, sizeof(struct FMMData)) ;   

  NbrElmsInSupport = 0 ;  
  NbrGroupsInSupport = List_Nbr(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList ); 

  iFMMGroup = 0 ;
  iFMMParentGroup = 0 ;
  
  for (iGroup = 0 ; iGroup < NbrGroupsInSupport ; iGroup++){
    List_Read(((struct Group *) List_Pointer(Problem_S.Group, InSupport))->InitialList, iGroup, &Current.Region) ;    
    Get_ValueOfExpressionByIndex(Current.FMM.DivXYZIndex, NULL, 0., 0., 0., &DivXYZ) ;
    NbrDivX0 = (int)DivXYZ.Val[0] ;
    NbrDivY0 = (int)DivXYZ.Val[1] ;
    NbrDivZ0 = (int)DivXYZ.Val[2] ;

    NbrDivX1 = (int)DivXYZ.Val[MAX_DIM] ;
    NbrDivY1 = (int)DivXYZ.Val[MAX_DIM+1] ;
    NbrDivZ1 = (int)DivXYZ.Val[MAX_DIM+2] ;

    Msg(INFO,"Div. Parents -> Region = %d DivX0 = %d  DivY0 = %d  DivZ0 = %d", Current.Region, NbrDivX0, NbrDivY0, NbrDivZ0) ; 
    Msg(INFO,"Div. Children -> DivX1 = %d  DivY1 = %d  DivZ1 = %d", NbrDivX1, NbrDivY1, NbrDivZ1) ; 
    
    NbrElmsInRegion = 0 ; 
    Xmax = Ymax = Zmax = Xmin = Ymin = Zmin = 0. ;

    for(iElm = 0 ; iElm < NbrElms ; iElm++){
      if((Elements_P0+iElm)->Region == Current.Region){
	NbrElmsInRegion++ ; 
	Geo_SetXYZMinMaxInElement((Elements_P0+iElm)->Num, &Xmin, &Ymin, &Zmin, &Xmax, &Ymax, &Zmax ) ;      
      }      
    }

    NbrElmsInSupport += NbrElmsInRegion ; 

    NbrElmsDiv = (int ***)Malloc(NbrDivX*sizeof(int**));
    for (j = 0 ; j < NbrDivX ; j++){
      NbrElmsDiv[j] = (int **)Calloc(NbrDivY, sizeof(int*));
      for (k = 0 ; k < NbrDivY ; k++)
	NbrElmsDiv[j][k] = (int *)Calloc(NbrDivZ, sizeof(int));
    }
    
    ElmsDiv = (int ****)Malloc(NbrDivX*sizeof(int***));
    for (j = 0 ; j < NbrDivX ; j++){
      ElmsDiv[j] = (int ***)Malloc(NbrDivY*sizeof(int**)); 
      for (k = 0 ; k < NbrDivY ; k++){
	ElmsDiv[j][k] = (int **)Malloc(NbrDivZ*sizeof(int*));
	for (l = 0 ; l < NbrDivZ ; l++)
	  ElmsDiv[j][k][l] = (int *)Malloc(sizeof(int));
      }
    }


    NbrGrsDiv = (int ***)Malloc(NbrDivX0*sizeof(int**));
    for (j = 0 ; j < NbrDivX0 ; j++){
      NbrGrsDiv[j] = (int **)Calloc(NbrDivY0, sizeof(int*));
      for (k = 0 ; k < NbrDivY0 ; k++)
	NbrGrsDiv[j][k] = (int *)Calloc(NbrDivZ0, sizeof(int));
    }


    GrsDiv = (int ****)Malloc(NbrDivX0*sizeof(int***));
    for (j = 0 ; j < NbrDivX0 ; j++){
      GrsDiv[j] = (int ***)Malloc(NbrDivY0*sizeof(int**)); 
      for (k = 0 ; k < NbrDivY0 ; k++){
	GrsDiv[j][k] = (int **)Malloc(NbrDivZ0*sizeof(int*));
	for (l = 0 ; l < NbrDivZ0 ; l++)
	  GrsDiv[j][k][l] = (int *)Malloc(sizeof(int));
      }
    }


    NbrDivX = NbrDivX0 * NbrDivX1 ;
    NbrDivY = NbrDivY0 * NbrDivY1 ;
    NbrDivZ = NbrDivZ0 * NbrDivZ1 ;

    for (iElm = 0 ; iElm < NbrElms ; iElm++){ 
      FoundHome = 0;
      
      for (j = 0 ; j < NbrDivX ; j++){
	for (k = 0 ; k < NbrDivY ; k++){
	  for (l = 0 ; l < NbrDivZ ; l++){
	    x1 = -tol+Xmin +(double)j*(Xmax-Xmin)/(double)NbrDivX;
	    x2 =  tol+Xmin +(double)(j+1)*(Xmax-Xmin)/(double)NbrDivX;
	    y1 = -tol+Ymin +(double)k*(Ymax-Ymin)/(double)NbrDivY;
	    y2 =  tol+Ymin +(double)(k+1)*(Ymax-Ymin)/(double)NbrDivY;
	    z1 = -tol+Zmin +(double)l*(Zmax-Zmin)/(double)NbrDivZ;
	    z2 =  tol+Zmin +(double)(l+1)*(Zmax-Zmin)/(double)NbrDivZ;

	    if( Current.Region == (Elements_P0+iElm)->Region ){
	      Geo_SetCentroidCoordinates((Elements_P0+iElm)->Num, Centroid);
	      
	      if (x1 <= Centroid[0] && Centroid[0] <= x2 &&
		  y1 <= Centroid[1] && Centroid[1] <= y2 &&
		  z1 <= Centroid[2] && Centroid[2] <= z2 ){
		if(!NbrElmsDiv[j][k][l]){
		  NbrElmsDiv[j][k][l] = 1;
		  ElmsDiv[j][k][l][0] = iElm ;
		} else {
		  NbrElmsDiv[j][k][l] += 1;
		  ElmsDiv[j][k][l] = (int *)realloc(ElmsDiv[j][k][l], NbrElmsDiv[j][k][l]*sizeof(int));
		  ElmsDiv[j][k][l][NbrElmsDiv[j][k][l]-1] = iElm;	    
		}
		FoundHome = 1; break;
	      }
	    }/* if Region */
	    if(FoundHome)break;
	  }
	  if(FoundHome)break;
	}
	if(FoundHome)break;
      }
      if(!FoundHome && (Current.Region == (Elements_P0+iElm)->Region)) 
	Msg(WARNING,"Element %d has not been assigned to a FMMGroup", iElm);
    }/* Elms NbrElms */ 


    NbrFMMGroupsInRegion = 0;
    for (j = 0 ; j < NbrDivX ; j++)
      for (k = 0 ; k < NbrDivY ; k++)
	for (l = 0 ; l < NbrDivZ ; l++)
	  if(NbrElmsDiv[j][k][l]){
	    Geo_InitFMMData(&FMMData_S) ;
	    FMMData_S.Element = List_Create(NbrElmsDiv[j][k][l], 1, sizeof(int)) ;
	    for(m = 0 ; m < NbrElmsDiv[j][k][l] ; m++){
	      List_Add(FMMData_S.Element, &(Elements_P0 + ElmsDiv[j][k][l][m])->Num) ;
	      Geo_SetCentroidCoordinates((Elements_P0+ElmsDiv[j][k][l][m])->Num, Centroid) ;

	      FMMData_S.Xgc += Centroid[0] ;
	      FMMData_S.Ygc += Centroid[1] ;
	      FMMData_S.Zgc += Centroid[2] ;
	      
	      (Elements_P0 + ElmsDiv[j][k][l][m])->FMMGroup = iFMMGroup ;
	    }	
	    
	    FMMData_S.Group = iFMMGroup;  
	    FMMData_S.Xgc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Ygc /= NbrElmsDiv[j][k][l];
	    FMMData_S.Zgc /= NbrElmsDiv[j][k][l];

	    List_Add(FMMGroup_L, &FMMData_S) ;
	    
	    iFMMGroup +=1 ;
	    NbrFMMGroupsInRegion += 1 ;
	  }/* if( NbrElmsDiv[j][k][l] ) */ 
  
    Free(NbrElmsDiv);
    Free(ElmsDiv); 

    Msg(INFO,"SupportIndex = %d  Region = %d  NbrElmsInRegion = %d  NbrFMMGroupsInRegion = %d",
	InSupport, Current.Region, NbrElmsInRegion, NbrFMMGroupsInRegion) ;
  }/* iGroup NbrGroupInSuppport*/
  
  Msg(INFO,"SupportIndex = %d  NbrElmsInSupport = %d  NbrFMMGroupsInSupport = %d",
      InSupport, NbrElmsInSupport, iFMMGroup) ; 

  for (iFMMGroup = 0 ; iFMMGroup < List_Nbr(FMMGroup_L) ; iFMMGroup++){       
    FoundHome = 0;
    List_Read(FMMGroup_L, iFMMGroup, &FMMData_S); 
    for (j = 0 ; j < NbrDivX0 ; j++){
      for (k = 0 ; k < NbrDivY0 ; k++){
	for (l = 0 ; l < NbrDivZ0 ; l++){
	  x1 = -tol+Xmin +(double)j*(Xmax-Xmin)/(double)NbrDivX0;
	  x2 =  tol+Xmin +(double)(j+1)*(Xmax-Xmin)/(double)NbrDivX0;
	  y1 = -tol+Ymin +(double)k*(Ymax-Ymin)/(double)NbrDivY0;
	  y2 =  tol+Ymin +(double)(k+1)*(Ymax-Ymin)/(double)NbrDivY0;
	  z1 = -tol+Zmin +(double)l*(Zmax-Zmin)/(double)NbrDivZ0;
	  z2 =  tol+Zmin +(double)(l+1)*(Zmax-Zmin)/(double)NbrDivZ0;
	  
	  if (x1 <=  FMMData_S.Xgc && FMMData_S.Xgc <= x2 &&
	      y1 <=  FMMData_S.Ygc && FMMData_S.Ygc <= y2 &&
	      z1 <=  FMMData_S.Zgc && FMMData_S.Zgc <= z2 ){
	    if(!NbrGrsDiv[j][k][l]){
	      NbrGrsDiv[j][k][l] = 1;
	      GrsDiv[j][k][l][0] = iFMMGroup ;
	    }
	    else {
	      NbrGrsDiv[j][k][l] += 1;
	      GrsDiv[j][k][l] = (int *)realloc(GrsDiv[j][k][l], NbrGrsDiv[j][k][l]*sizeof(int));
	      GrsDiv[j][k][l][NbrGrsDiv[j][k][l]-1] = iFMMGroup;
	    }
	    FoundHome = 1; break;
	  }	    
	  if(FoundHome)break;
	}
	if(FoundHome)break;
      }
      if(FoundHome)break;
    }
  } /* Grs in List_Nbr(FMMGroup_L) */ 
  
  for (j = 0 ; j < NbrDivX0 ; j++)
    for (k = 0 ; k < NbrDivY0 ; k++)
      for (l = 0 ; l < NbrDivZ0 ; l++)
	if(NbrGrsDiv[j][k][l]){
	  Geo_InitFMMData(&FMMParentData_S) ;
	  FMMParentData_S.Element = List_Create(NbrGrsDiv[j][k][l], 1, sizeof(struct FMMData)) ;
	  for(m = 0 ; m < NbrGrsDiv[j][k][l] ; m++){
	    List_Read(FMMGroup_L, GrsDiv[j][k][l][m], &FMMData_S); 
	    List_Add(FMMParentData_S.Element, &FMMData_S) ;
	    
	    FMMParentData_S.Xgc +=  FMMData_S.Xgc ;
	    FMMParentData_S.Ygc +=  FMMData_S.Ygc ;
	    FMMParentData_S.Zgc +=  FMMData_S.Zgc ;   
	  }	
	  
	  FMMParentData_S.Group = iFMMParentGroup;  
	  FMMParentData_S.Xgc /= NbrGrsDiv[j][k][l];
	  FMMParentData_S.Ygc /= NbrGrsDiv[j][k][l];
	  FMMParentData_S.Zgc /= NbrGrsDiv[j][k][l];
	  
	  List_Add(FMMParentGroup_L, &FMMParentData_S) ;
	  
	  iFMMParentGroup +=1 ;
	}/* if( NbrGrsDiv[j][k][l] ) */
  
  
  Free(NbrGrsDiv);
  Free(GrsDiv);
  
  List_Reset(FMMGroup_L);
  
  FMMGroup_S.InIndex = InSupport ;
  FMMGroup_S.List = FMMParentGroup_L ;
  List_Add(Problem_S.FMMGroup, &FMMGroup_S ); 
  
  GetDP_End ;
}

#endif


void Get_InFMMGroupList( int Index_Formulation, struct GeoData *GeoData_P ){

  int  i_EquTerm, Nbr_EquationTerm, i_Observation, i_Source, FlagError, DIM ;
  double k0 = -1. ;
  struct Formulation    * Formulation_P ;
  struct EquationTerm   * EquationTerm_P0, * EquationTerm_P ;
  struct DefineQuantity * DefineQuantityDof_P, * DefineQuantity_P0 ;
  struct FMMmat         * FMMmat_P, FMMmat_S ;
  List_T *InIndex ;

  GetDP_Begin("Get_InFMMGroupList") ;

  Formulation_P = (struct Formulation*)List_Pointer(Problem_S.Formulation, Index_Formulation) ;
  EquationTerm_P0 = (struct EquationTerm*)List_Pointer(Formulation_P->Equation, 0) ;
  Nbr_EquationTerm = List_Nbr(Formulation_P->Equation) ;
  DefineQuantity_P0 = (struct DefineQuantity*)List_Pointer(Formulation_P->DefineQuantity, 0) ;

  InIndex = List_Create( 1, 1, sizeof(int)) ;

  Problem_S.FMMGroup = List_Create (1, 1, sizeof(struct FMMGroup)) ;

  Current.DofData->FMM_Matrix = NULL ;

  /* Current.FMM.Dimension = DIM = (int)((GeoData_P)->Dimension) ; */
  
  for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
    EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
    EquationTerm_P->Case.LocalTerm.FMMObservation =
      EquationTerm_P->Case.LocalTerm.FMMSource = EquationTerm_P->Case.LocalTerm.iFMMEqu = -1 ;

    if ( EquationTerm_P->Type == GALERKIN && EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex > 1 ) {
      DefineQuantityDof_P = DefineQuantity_P0 + EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof ;
      if( DefineQuantityDof_P->Type==INTEGRALQUANTITY){

	if(DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.NbrParameters == 2)
	  k0 = DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Para[1] ;/* Helmoltz Case */
	  DIM = (int)DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Para[0] ;
	  Current.FMM.Dimension = DIM ;

	if( List_PQuery(InIndex, &EquationTerm_P->Case.LocalTerm.InIndex, fcmp_int) == NULL ){
	  List_Add( InIndex, &EquationTerm_P->Case.LocalTerm.InIndex ) ;
	  Geo_CreateFMMGroup( EquationTerm_P->Case.LocalTerm.InIndex, GeoData_P, k0) ;
	}

	if( List_PQuery(InIndex, &DefineQuantityDof_P->IntegralQuantity.InIndex, fcmp_int) == NULL ){
	  List_Add( InIndex, &DefineQuantityDof_P->IntegralQuantity.InIndex ) ;
	  Geo_CreateFMMGroup( DefineQuantityDof_P->IntegralQuantity.InIndex, GeoData_P, k0 ) ;
	}

	

	if (Current.DofData->FMM_Matrix == NULL)
	  Current.DofData->FMM_Matrix = List_Create(1, 1, sizeof(struct FMMmat)) ;
	
	i_Source = EquationTerm_P->Case.LocalTerm.FMMSource = 
	  List_ISearch(InIndex, &DefineQuantityDof_P->IntegralQuantity.InIndex, fcmp_int) ;
	i_Observation = EquationTerm_P->Case.LocalTerm.FMMObservation = 
	  List_ISearch(InIndex, &EquationTerm_P->Case.LocalTerm.InIndex, fcmp_int) ;

	InitFMMmatrix(i_Source, i_Observation, &FMMmat_S) ;
	List_Add( Current.DofData->FMM_Matrix, &FMMmat_S) ;
 	EquationTerm_P->Case.LocalTerm.iFMMEqu = List_Nbr(Current.DofData->FMM_Matrix)-1 ;
	
	Get_GroupNeighbours( EquationTerm_P->Case.LocalTerm.iFMMEqu) ;

	FMMmat_P = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, EquationTerm_P->Case.LocalTerm.iFMMEqu) ;

	if (DIM == _2D)
	  Get_FunctionForFunction( FMMProd_Function2D, DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Fct,
				   &FlagError, &(FMMmat_P->FctProd) ) ;
	else
	  Get_FunctionForFunction( FMMProd_Function3D, DefineQuantityDof_P->IntegralQuantity.FunctionForFMM.Fct,
				   &FlagError, &(FMMmat_P->FctProd) ) ;

	FMMmat_P->GFx = &(DefineQuantityDof_P->IntegralQuantity.FunctionForFMM) ;

	
      }/* Quantity */
    }/* If Galerkin */
  }/* Equation_Term */

    
  Init_CurrentFMMData ( k0 ) ;
  List_Delete(InIndex);

  GetDP_End ; 

}