fmm_caldtamatrices.c

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

#define RCSID "$Id: FMM_CalDTAmatrices.c,v 1.11 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 "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "CurrentData.h"
#include "Get_Geometry.h"
#include "Get_DofOfElement.h"
#include "DofData.h"
#include "Data_DefineE.h"
#include "Cal_Value.h"
#include "GF_Function.h"
#include "Tools.h"
#include "Data_FMM.h"
#include "Cal_FMMAnalyticalIntegral.h"
#include "F_FMM_DTA.h"

/* ------------------------------------------------------------------------- */ 
/* TRANSLATION in the spectral domain                                        */ 
/* ------------------------------------------------------------------------- */ 

void GF_FMMTranslationValue ( ){
   
  int FMMObs, FMMSrc, NbrGroupSrc, NbrFG, *FG ;
  int i, j, Nd, iEqu, NbrFMMEqu, jEqu, N, NbrHar ;
  double ***T ;

  struct Value TV[NBR_MAX_DIR] ;
  struct FMMData *FMMDataObs_P0,  *FMMDataSrc_P0 ;
  struct FMMGroup  FMMGroup_S ;  
  struct FMMmat *FMMmat_P0 ;
  struct Function *GFx ;
  List_T *FG_L ;


  GetDP_Begin("GF_FMMTranslationValue");

  Msg(INFO, "Creation Translation matrix (with GF_Functions)");

  Nd = Current.FMM.NbrDir ;
 
  /* The dimensions of the matrix depend on the case we are dealing with.
     Helmholtz -> Nd, Laplace2D -> Nd * Nd, Laplace3D -> (2*Nd-1)*(2*Nd+1) */

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

  for(iEqu = 0 ; iEqu < NbrFMMEqu ; iEqu ++ ){
    FMMSrc = (FMMmat_P0+iEqu)->Src ;
    FMMObs = (FMMmat_P0+iEqu)->Obs ;
    GFx = (FMMmat_P0+iEqu)->GFx ;

    N = (GFx->NbrParameters==2 || GFx->NbrParameters == 3 ) ? Nd : 
      (((int)GFx->Para[0] == _2D) ? Nd*Nd : (2*Nd-1)*(2*Nd+1)) ; 

    jEqu = 0 ; while( jEqu<iEqu && ((FMMmat_P0+jEqu)->Src != FMMSrc ||
				    (FMMmat_P0+jEqu)->Obs != FMMObs)) jEqu ++ ;
    if (jEqu < iEqu && GFx->NbrParameters == 1 && GFx->Para[0] != _2D  )
      (FMMmat_P0 + iEqu)->T = (FMMmat_P0 + jEqu)->T ;
    else { 
      /* If the Source and Observation supports are the same for two Galerkin Terms, 
	 the translation matrix is common */ 
      List_Read(Problem_S.FMMGroup, FMMSrc, &FMMGroup_S ) ;
      NbrGroupSrc = List_Nbr(FMMGroup_S.List ) ;
      FMMDataSrc_P0 = FMMDataObs_P0 =(struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;
      
      if ( FMMSrc != FMMObs ){
	List_Read(Problem_S.FMMGroup, FMMObs, &FMMGroup_S ) ;
	FMMDataObs_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List,0) ;   
      }
 
      T = (double***)Malloc(NbrGroupSrc*sizeof(double**)) ;
      for(i = 0 ; i < NbrGroupSrc ; i++){
	List_Read((FMMmat_P0+iEqu)->FG_L, i, &FG_L) ;	  
	NbrFG = List_Nbr(FG_L) ;
	T[i] = (double**)Malloc(NbrFG*sizeof(double*)) ;
	for (j = 0 ; j < NbrFG ; j++)
	    T[i][j] = (double*)Malloc(2*N*sizeof(double)) ;	
      }
      
      Current.FMM.Flag_GF = FMM_TRANSLATION ;

      for(i = 0 ; i < NbrGroupSrc ; i++){ /* i Origin */
	List_Read((FMMmat_P0+iEqu)->FG_L, i, &FG_L) ;	  
	NbrFG = List_Nbr(FG_L) ;
	if(NbrFG){
	  FG = (int*)(FG_L->array) ;
	  Current.xs = (FMMDataSrc_P0+i)->Xgc ;
	  Current.ys = (FMMDataSrc_P0+i)->Ygc ;
	  Current.zs = (FMMDataSrc_P0+i)->Zgc ;
	  Current.FMM.Rsrc = (FMMDataSrc_P0+i)->Rmax ;

	  for (j = 0 ; j < NbrFG ; j++){ /* j Destination */
	    Current.x = (FMMDataObs_P0 + FG[j])->Xgc ;
	    Current.y = (FMMDataObs_P0 + FG[j])->Ygc ;
	    Current.z = (FMMDataObs_P0 + FG[j])->Zgc ;
	    Current.FMM.Robs = (FMMDataObs_P0 + FG[j])->Rmax ;
	    ((void(*)(struct Function*, struct Value*, struct Value*))GFx->Fct)(GFx, TV, TV) ;
	    Cal_ValueArray2DoubleArray(TV, T[i][j], N) ;
	  }
	}
       	else T[i] = NULL ;
      }  
      (FMMmat_P0+iEqu)->T = T ;
    }
  }
  Current.NbrHar = NbrHar ;
  Current.FMM.Flag_GF = FMM_DIRECT ;

  GetDP_End ;
}



void GF_FMMTranslationValueAdaptive( ){
   
  int FMMObs, FMMSrc, NbrGroupSrc, NbrFG, *FG ;
  int i, j, iEqu, NbrFMMEqu, jEqu, N, *Nd_A, NbrDirMAX, NbrHar ;
  double ***T ;

  struct Value TV[NBR_MAX_DIR] ;
  struct FMMData *FMMDataObs_P0,  *FMMDataSrc_P0 ;
  struct FMMGroup  FMMGroup_S ;  
  struct FMMmat *FMMmat_P0 ;
  struct Function *GFx ;
  List_T *FG_L, *Nd_L ;


  GetDP_Begin("GF_FMMTranslationValueAdaptive");
  Msg(RESOURCES, "Before translation ");
  Msg(INFO, "Creation Translation matrix (with GF_Functions & Adaptive)");
 
  /* ADAPTIVE number of directions: Truncation parameter is a function of
     the precision and the ratio D/R, where D is the distance between SRC
     group and OBS group and R is the maximum radius of the circles/spheres
     enclosing the groups SRC and OBS. 
     The Aggregation and Disaggregation matrices are built for the extreme case,
     the suitable number of elements will be taken into account during the
     iterative process.
     The dimensions of the matrix depend on the case we are dealing with.
     Helmholtz -> Nd, Laplace2D -> Nd * Nd, Laplace3D -> (2*Nd-1)*(2*Nd+1) */

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


  for(iEqu = 0 ; iEqu < NbrFMMEqu ; iEqu ++ ){
    FMMSrc = (FMMmat_P0+iEqu)->Src ;
    FMMObs = (FMMmat_P0+iEqu)->Obs ;
    GFx = (FMMmat_P0+iEqu)->GFx ;

    jEqu = 0 ; while( jEqu<iEqu && ((FMMmat_P0+jEqu)->Src != FMMSrc ||
				    (FMMmat_P0+jEqu)->Obs != FMMObs)) jEqu ++ ;
    if ( jEqu < iEqu  && GFx->NbrParameters == 1 && GFx->Para[0] != _2D ){
      (FMMmat_P0 + iEqu)->T = (FMMmat_P0 + jEqu)->T ;
    }
    /* If the Source and Observation supports are the same for two Galerkin Terms, 
       the translation matrix is common */ 
    else {
      List_Read(Problem_S.FMMGroup, FMMSrc, &FMMGroup_S ) ;
      NbrGroupSrc = List_Nbr(FMMGroup_S.List ) ;
      FMMDataSrc_P0 = FMMDataObs_P0 =(struct FMMData*)List_Pointer(FMMGroup_S.List, 0) ;
      
      if ( FMMSrc != FMMObs ){
	List_Read(Problem_S.FMMGroup, FMMObs, &FMMGroup_S ) ;
	FMMDataObs_P0 = (struct FMMData*)List_Pointer(FMMGroup_S.List,0) ;   
      }
 
      T = (double***)Malloc(NbrGroupSrc*sizeof(double**)) ;
      for(i = 0 ; i < NbrGroupSrc ; i++){
	List_Read((FMMmat_P0+iEqu)->FG_L, i, &FG_L) ;	  
	NbrFG = List_Nbr(FG_L) ;
	List_Read((FMMmat_P0+iEqu)->Nd_L, i, &Nd_L) ;
	Nd_A = (int*)(Nd_L->array);
	T[i] = (double**)Malloc(NbrFG*sizeof(double*)) ;
	
	for (j = 0 ; j < NbrFG ; j++){
	  N = (GFx->NbrParameters==2) ? Nd_A[j] :
	    (((int)GFx->Para[0] == _2D) ? Nd_A[j]*Nd_A[j] : (2*Nd_A[j]-1)*(2*Nd_A[j]+1) ); 
	  T[i][j] = (double*)Malloc(2*N*sizeof(double)) ;
	}	
      }
      
      
      Current.FMM.Flag_GF = FMM_TRANSLATION ;

      for(i = 0 ; i < NbrGroupSrc ; i++){ /* i Origin */
	List_Read((FMMmat_P0+iEqu)->FG_L, i, &FG_L) ;	  
	List_Read((FMMmat_P0+iEqu)->Nd_L, i, &Nd_L) ;	  
	NbrFG = List_Nbr(FG_L) ;
	if(NbrFG){
	  FG = (int*)(FG_L->array) ;
	  Nd_A = (int*)(Nd_L->array) ;
	  Current.xs = (FMMDataSrc_P0+i)->Xgc ;
	  Current.ys = (FMMDataSrc_P0+i)->Ygc ;
	  Current.zs = (FMMDataSrc_P0+i)->Zgc ;
	  Current.FMM.Rsrc = (FMMDataSrc_P0+i)->Rmax ;

	  for (j = 0 ; j < NbrFG ; j++){ /* j Destination */
	    Current.FMM.NbrDir = Nd_A[j] ;
	    N = (GFx->NbrParameters==2) ? Nd_A[j] :
	      (((int)GFx->Para[0] == _2D) ? Nd_A[j]*Nd_A[j] : (2*Nd_A[j]-1)*(2*Nd_A[j]+1) );
	    Current.x = (FMMDataObs_P0 + FG[j])->Xgc ;
	    Current.y = (FMMDataObs_P0 + FG[j])->Ygc ;
	    Current.z = (FMMDataObs_P0 + FG[j])->Zgc ;
	    Current.FMM.Robs = (FMMDataObs_P0 + FG[j])->Rmax ;
	    
	    ((void(*)(struct Function*, struct Value*, struct Value*)) GFx->Fct)(GFx, TV, TV) ;
	    
	    Cal_ValueArray2DoubleArray(TV, T[i][j], N) ;
	  }
	}
       	else T[i] = NULL ;
      } 
      (FMMmat_P0+iEqu)->T = T ;
    }
  }

  Current.FMM.NbrDir = NbrDirMAX ;
  Current.NbrHar = NbrHar ;
  Current.FMM.Flag_GF = FMM_DIRECT ;

  Msg(RESOURCES, "After translation");

  GetDP_End ;
}



#define CAST3V    void(*)(struct Value*, struct Value*, struct Value*)
#define CASTF2V   void(*)(struct Function*, struct Value*, struct Value*)

/*-------------------------------------------------------------------------------*/
/*           C a l _ F M M G a l e r k i n D i s a g g r e g a t i o n           */
/*-------------------------------------------------------------------------------*/

void  Cal_FMMGalerkinDisaggregation(struct EquationTerm     * EquationTerm_P0,
				    struct QuantityStorage  * QuantityStorage_P0) {
  
  struct EquationTerm       * EquationTerm_P ;
  struct FemLocalTermActive * FI ;
  struct QuantityStorage    * QuantityStorageEqu_P, * QuantityStorageDof_P;
  struct IntegrationCase    * IntegrationCase_P ;
  struct Quadrature         * Quadrature_P ;
  struct Value                vBFxEquV ;
  struct Value                TermFct ;
  struct Value                GFValue [NBR_MAX_DIR],  BFGFValue [NBR_MAX_DIR] ;
  struct Element              Element ;
  struct FMMData            * FMMData_P0, * FMMData_P ; 
  struct FMMGroup             FMMGroup_S ;  
  struct FMMmat             * FMMmat_P0 , * FMMmat_P ;
  double                   ** Disagg_M ;
  List_T                    * NumEqu_L ;

  int     Nbr_Equ,  FMMObs, NbrFMMEqu, i_FMMEqu, NbrDir, i_FMM ;
  int     Nbr_Group, i_Group, Nbr_ElmsInGroup, i_Element, *ElmtsGr ;
  int     i, j, Type_Dimension, Nbr_IntPoints, i_IntPoint ;
  int     Nbr_EquList, NbrHar ;
  
  double  weight, Factor ;
  double  vBFuEqu [MAX_DIM], vBFxEqu [MAX_DIM] ;


  void (*xFunctionBFEqu[NBR_MAX_BASISFUNCTIONS])
    (struct Element * Element, int NumEntity, 
     double u, double v, double w, double Value[] ) ;
  double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
  void (*Get_IntPoint)(int,int,double*,double*,double*,double*);

  List_T * FMMIntegrationCaseEqu_L ;
  int      FMMCriterionIndexEqu ;


  GetDP_Begin("Cal_FMMGalerkinDisaggregation");