gf_helmholtz.c

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

#define RCSID "$Id: GF_Helmholtz.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 <[e-mail]>.
 *
 * Contributor(s):
 *   Ruth Sabariego
 */

#include "GetDP.h" 
#include "Data_Active.h"
#include "Cal_Value.h"
#include "CurrentData.h"
#include "Numeric.h"
#include "Data_FMM.h"
#include "Legendre.h"
#include "SphBessel.h"


/* ------------------------------------------------------------------------ */
/*  Warning: the pointers A and V can be identical. You must                */
/*           use temporary variables in your computations: you can only     */
/*           affect to V at the very last time (when you're sure you will   */
/*           not use A any more).                                           */
/* ------------------------------------------------------------------------ */

#define F_ARG     struct Function * Fct, struct Value * A, struct Value * V

/* 
   Fct->Para[0] == dimension 
   Fct->Para[1] == wave number
*/



/* ------------------------------------------------------------------------ */
/*  G F _ H e l m h o l t z                                                 */
/* ------------------------------------------------------------------------ */

void GF_Helmholtz (F_ARG) {
 
  double  r, phi_r, phi, c, s, kr, k0r, bjn, byn, cr ;
  double xxs, yys, zzs, sbjn, sbyn, Pl ;
  double P[NBR_MAX_EXP], jsph[NBR_MAX_EXP], ysph[NBR_MAX_EXP] ;
  int ns, Nd, i_FMM, j_FMM ;

  GetDP_Begin("GF_Helmholtz");

  if(Current.NbrHar != 2)
    Msg(GERROR, "Wrong Number of Harmonics in 'GF_Helmholtz'");
  

  V->Type = SCALAR ;
  
  switch((int)Fct->Para[0]){      
  case _2D :
    switch(Current.FMM.Flag_GF){            
    case FMM_DIRECT :
      r = sqrt(SQU(Current.x-Current.xs)+
	       SQU(Current.y-Current.ys) ) ;
      if(!r) Msg(GERROR, "1/0 in 'GF_Helmholtz'") ;
      kr = Fct->Para[1]*r;

      V->Val[0]       = -y0(kr)/4 ; 
      V->Val[MAX_DIM] = -j0(kr)/4 ;
      break ;
      
    case FMM_AGGREGATION :
      Nd = Current.FMM.NbrDir ;
      
      r = sqrt( SQU(Current.xs-Current.FMM.Xgc)+ SQU(Current.ys-Current.FMM.Ygc)) ;
      phi_r = atan2(Current.ys- Current.FMM.Ygc, Current.xs-Current.FMM.Xgc) ;     
      kr = Fct->Para[1]*r;
      
      if (r==0) 
	for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	  V[i_FMM].Type = SCALAR ;
	  V[i_FMM].Val[0] = 1. ;
	  V[i_FMM].Val[MAX_DIM] = 0. ;
	}
      else
	for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	  phi = Current.FMM.Phi[i_FMM] ;
	  V[i_FMM].Type = SCALAR ;
	  V[i_FMM].Val[0]       =  cos(kr*sin(phi+phi_r)) ;
	  V[i_FMM].Val[MAX_DIM] = -sin(kr*sin(phi+phi_r)) ;	  
	}
      break;

    case FMM_DISAGGREGATION :
      Nd =  Current.FMM.NbrDir ;

      r = sqrt( SQU(Current.x - Current.FMM.Xgc)+ SQU(Current.y - Current.FMM.Ygc)) ;
      phi_r = atan2(Current.y - Current.FMM.Ygc, Current.x - Current.FMM.Xgc) ;
      kr = Fct->Para[1]*r;

      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	V[i_FMM].Type = SCALAR ;
	phi = Current.FMM.Phi[i_FMM];
	V[i_FMM].Val[0]       = cos(kr*sin(phi+phi_r)) ;
	V[i_FMM].Val[MAX_DIM] = sin(kr*sin(phi+phi_r)) ;
      }
      break;

    case FMM_TRANSLATION :
      Nd =  Current.FMM.NbrDir ;
      ns = (Nd-1)/2 ;
      cr = 1./4./Nd ; 

      kr = Fct->Para[1] * sqrt( SQU(Current.x - Current.xs)+ SQU(Current.y - Current.ys)) ;
      phi_r = atan2(Current.y - Current.ys, Current.x - Current.xs) ;
      
      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){
	V[i_FMM].Type = SCALAR ;
	phi = Current.FMM.Phi[i_FMM];
	V[i_FMM].Val[0]       = - cr*yn(0,kr);
	V[i_FMM].Val[MAX_DIM] = - cr*jn(0,kr);

	for (j_FMM = 1 ; j_FMM <= ns ;  j_FMM++){ 
	  c = cos(j_FMM*(phi+phi_r-PI));
	  s = sin(j_FMM*(phi+phi_r-PI));
	  bjn = jn(j_FMM,kr);
	  byn = yn(j_FMM,kr);
	  V[i_FMM].Val[0]       += -2.*cr* (j_FMM%2 ?  bjn*s : byn*c) ;
	  V[i_FMM].Val[MAX_DIM] += -2.*cr* (j_FMM%2 ? -byn*s : bjn*c) ; 
	}
      }
      break;      

    default :
      Msg(GERROR, "Case 2D: Bad Flag_GF 'GF_Helmholtz' (%d)", Current.FMM.Flag_GF);  
      break;
    }
    break;

  case _3D :
    switch(Current.FMM.Flag_GF){          
    case FMM_DIRECT :
      r = sqrt(SQU(Current.x-Current.xs)+
	       SQU(Current.y-Current.ys)+
	       SQU(Current.z-Current.zs)) ;
      if(!r) Msg(GERROR, "1/0 in 'GF_Helmholtz'") ;

      kr = Fct->Para[1]*r;     
      V->Val[0]       =  ONE_OVER_FOUR_PI * cos(kr) / r ; 
      V->Val[MAX_DIM] = -ONE_OVER_FOUR_PI * sin(kr) / r ;
      break ;

    case FMM_AGGREGATION : /* 3D GF_Helmholtz */
      Nd = Current.FMM.NbrDir ;
      xxs = Current.xs-Current.FMM.Xgc ;
      yys = Current.ys-Current.FMM.Ygc ;
      zzs = Current.zs-Current.FMM.Zgc ;
      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	kr = Fct->Para[1] * (xxs*Current.FMM.Kdir[i_FMM][0] +
			     yys*Current.FMM.Kdir[i_FMM][1] +
			     zzs*Current.FMM.Kdir[i_FMM][2] ) ;	
	V[i_FMM].Type = SCALAR ;
	V[i_FMM].Val[0]       =  cos(kr) ;
	V[i_FMM].Val[MAX_DIM] = -sin(kr) ;
      }
      break;

    case FMM_DISAGGREGATION : /* 3D GF_Helmholtz*/
      Nd = Current.FMM.NbrDir ;
      xxs = Current.x-Current.FMM.Xgc ;
      yys = Current.y-Current.FMM.Ygc ;
      zzs = Current.z-Current.FMM.Zgc ;

      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	kr = Fct->Para[1] * (xxs*Current.FMM.Kdir[i_FMM][0] +
			     yys*Current.FMM.Kdir[i_FMM][1] +
			     zzs*Current.FMM.Kdir[i_FMM][2] ) ;	
	V[i_FMM].Type = SCALAR ;
	V[i_FMM].Val[0]       = cos(kr) ;
	V[i_FMM].Val[MAX_DIM] = sin(kr) ;	  
      }       
      break;

    case FMM_TRANSLATION : /* 3D GF_Helmholtz*/
      Nd =  Current.FMM.NbrDir ;
      ns =  (int)sqrt((double)Nd/2) ;
      cr = - Fct->Para[1] * ONE_OVER_FOUR_PI * ns/2/Nd ; 

      xxs = Current.xs - Current.x ;
      yys = Current.ys - Current.y ;
      zzs = Current.zs - Current.z ;
      r =  sqrt( SQU(xxs) + SQU(yys) + SQU(zzs) ) ; 
      k0r = Fct->Para[1] * r ;     

      Spherical_j_nArray(0, k0r, ns+1, jsph) ;
      Spherical_y_nArray(0, k0r, ns+1, ysph) ;
      
      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	V[i_FMM].Type = SCALAR ;
	V[i_FMM].Val[0] = V[i_FMM].Val[MAX_DIM] = 0. ;
	kr = (xxs*Current.FMM.Kdir[i_FMM][0] + 
	      yys*Current.FMM.Kdir[i_FMM][1] +
	      zzs*Current.FMM.Kdir[i_FMM][2])/r ; 

	LegendreRecursive( ns, 0, kr, P);
    	
	for (j_FMM = 0 ; j_FMM <= ns ; j_FMM++){ 
	  s = ((j_FMM/2)%2)? -1 : 1 ;	 
	  sbjn = jsph[j_FMM] ;
	  sbyn = ysph[j_FMM] ;
	  Pl = (2*j_FMM+1)* P[j_FMM] ;

	  V[i_FMM].Val[0]       += s * Pl * ((j_FMM%2) ?   sbjn : sbyn ) ;
	  V[i_FMM].Val[MAX_DIM] += s * Pl * ((j_FMM%2) ? (-sbyn): sbjn ) ;
	}

	V[i_FMM].Val[0]       *= cr*Current.FMM.Weight[i_FMM];
	V[i_FMM].Val[MAX_DIM] *= cr*Current.FMM.Weight[i_FMM];
      }

      break;      
    default :
      Msg(GERROR, "Case 3D: Bad Flag_GF 'GF_Helmholtz' (%d)", Current.FMM.Flag_GF);  
      break;

    }
    break;

  default :
    Msg(GERROR, "Bad Parameter for 'GF_Helmholtz' (%d)", (int)Fct->Para[0]);
    break;
  }
     
  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  G F _ H e l m h o l t z T h i n W i r e                                 */
/* ------------------------------------------------------------------------ */



void GF_HelmholtzThinWire (F_ARG) {
 
  double  a , r, phi_r, phi, c, s, kr, k0r, bjn, byn, cr ;
  double xxs, yys, zzs, sbjn, sbyn, Pl ;
  double P[NBR_MAX_EXP], jsph[NBR_MAX_EXP], ysph[NBR_MAX_EXP] ;
  int ns, Nd, i_FMM, j_FMM ;

  GetDP_Begin("GF_HelmholtzThinWire");

  if(Current.NbrHar != 2)
    Msg(GERROR, "Wrong Number of Harmonics in 'GF_HelmholtzThinWire'");
  

  V->Type = SCALAR ;
  
  switch((int)Fct->Para[0]){      
  case _2D :
    switch(Current.FMM.Flag_GF){            
    case FMM_DIRECT :
      a =  Fct->Para[2] ;
      r = sqrt(SQU(Current.x-Current.xs)+
	       SQU(Current.y-Current.ys)+SQU(a)) ;
      if(!r) Msg(GERROR, "1/0 in 'GF_HelmholtzThinWire'") ;
      kr = Fct->Para[1]*r;

      V->Val[0]       = -y0(kr)/4 ; 
      V->Val[MAX_DIM] = -j0(kr)/4 ;
      break ;
      
    case FMM_AGGREGATION :
      Nd = Current.FMM.NbrDir ;
      
      r = sqrt( SQU(Current.xs-Current.FMM.Xgc)+ SQU(Current.ys-Current.FMM.Ygc)) ;
      phi_r = atan2(Current.ys- Current.FMM.Ygc, Current.xs-Current.FMM.Xgc) ;     
      kr = Fct->Para[1]*r;
      
      if (r==0) 
	for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	  V[i_FMM].Type = SCALAR ;
	  V[i_FMM].Val[0] = 1. ;
	  V[i_FMM].Val[MAX_DIM] = 0. ;
	}
      else
	for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	  phi = Current.FMM.Phi[i_FMM] ;
	  V[i_FMM].Type = SCALAR ;
	  V[i_FMM].Val[0]       =  cos(kr*sin(phi+phi_r)) ;
	  V[i_FMM].Val[MAX_DIM] = -sin(kr*sin(phi+phi_r)) ;	  
	}
      break;

    case FMM_DISAGGREGATION :
      Nd =  Current.FMM.NbrDir ;

      r = sqrt( SQU(Current.x - Current.FMM.Xgc)+ SQU(Current.y - Current.FMM.Ygc)) ;
      phi_r = atan2(Current.y - Current.FMM.Ygc, Current.x - Current.FMM.Xgc) ;
      kr = Fct->Para[1]*r;

      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	V[i_FMM].Type = SCALAR ;
	phi = Current.FMM.Phi[i_FMM];
	V[i_FMM].Val[0]       = cos(kr*sin(phi+phi_r)) ;
	V[i_FMM].Val[MAX_DIM] = sin(kr*sin(phi+phi_r)) ;
      }
      break;

    case FMM_TRANSLATION :
      Nd =  Current.FMM.NbrDir ;
      ns = (Nd-1)/2 ;
      cr = 1./4./Nd ; 

      kr = Fct->Para[1] * sqrt( SQU(Current.x - Current.xs)+ SQU(Current.y - Current.ys)) ;
      phi_r = atan2(Current.y - Current.ys, Current.x - Current.xs) ;
      
      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){
	V[i_FMM].Type = SCALAR ;
	phi = Current.FMM.Phi[i_FMM];
	V[i_FMM].Val[0]       = - cr*yn(0,kr);
	V[i_FMM].Val[MAX_DIM] = - cr*jn(0,kr);

	for (j_FMM = 1 ; j_FMM <= ns ;  j_FMM++){ 
	  c = cos(j_FMM*(phi+phi_r-PI));
	  s = sin(j_FMM*(phi+phi_r-PI));
	  bjn = jn(j_FMM,kr);
	  byn = yn(j_FMM,kr);
	  V[i_FMM].Val[0]       += -2.*cr* (j_FMM%2 ?  bjn*s : byn*c) ;
	  V[i_FMM].Val[MAX_DIM] += -2.*cr* (j_FMM%2 ? -byn*s : bjn*c) ; 
	}
      }
      break;      

    default :
      Msg(GERROR, "Case 2D: Bad Flag_GF 'GF_HelmholtzThinWire' (%d)", Current.FMM.Flag_GF);  
      break;
    }
    break;

  case _3D :
    switch(Current.FMM.Flag_GF){          
    case FMM_DIRECT :
      a =  Fct->Para[2] ;

      r = sqrt(SQU(Current.x-Current.xs)+
	       SQU(Current.y-Current.ys)+
	       SQU(Current.z-Current.zs)+SQU(a)) ;
      if(!r) Msg(GERROR, "1/0 in 'GF_HelmholtzThinWire'") ;
      
      kr = Fct->Para[1]*r;     
      V->Val[0]       =  ONE_OVER_FOUR_PI * cos(kr) / r ; 
      V->Val[MAX_DIM] = -ONE_OVER_FOUR_PI * sin(kr) / r ;
      break ;

    case FMM_AGGREGATION : /* 3D GF_Helmholtz */
      Nd = Current.FMM.NbrDir ;
      xxs = Current.xs-Current.FMM.Xgc ;
      yys = Current.ys-Current.FMM.Ygc ;
      zzs = Current.zs-Current.FMM.Zgc ;
      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	kr = Fct->Para[1] * (xxs*Current.FMM.Kdir[i_FMM][0] +
			     yys*Current.FMM.Kdir[i_FMM][1] +
			     zzs*Current.FMM.Kdir[i_FMM][2] ) ;	
	V[i_FMM].Type = SCALAR ;
	V[i_FMM].Val[0]       =  cos(kr) ;
	V[i_FMM].Val[MAX_DIM] = -sin(kr) ;
      }
      break;

    case FMM_DISAGGREGATION : /* 3D GF_Helmholtz*/
      Nd = Current.FMM.NbrDir ;
      xxs = Current.x-Current.FMM.Xgc ;
      yys = Current.y-Current.FMM.Ygc ;
      zzs = Current.z-Current.FMM.Zgc ;

      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	kr = Fct->Para[1] * (xxs*Current.FMM.Kdir[i_FMM][0] +
			     yys*Current.FMM.Kdir[i_FMM][1] +
			     zzs*Current.FMM.Kdir[i_FMM][2] ) ;	
	V[i_FMM].Type = SCALAR ;
	V[i_FMM].Val[0]       = cos(kr) ;
	V[i_FMM].Val[MAX_DIM] = sin(kr) ;	  
      }       
      break;

    case FMM_TRANSLATION : /* 3D GF_Helmholtz*/
      Nd =  Current.FMM.NbrDir ;
      ns =  (int)sqrt((double)Nd/2) ;
      cr = - Fct->Para[1] * ONE_OVER_FOUR_PI * ns/2/Nd ; 

      xxs = Current.xs - Current.x ;
      yys = Current.ys - Current.y ;
      zzs = Current.zs - Current.z ;
      r =  sqrt( SQU(xxs) + SQU(yys) + SQU(zzs) ) ; 
      k0r = Fct->Para[1] * r ;     

      Spherical_j_nArray(0, k0r, ns+1, jsph) ;
      Spherical_y_nArray(0, k0r, ns+1, ysph) ;
      
      for (i_FMM = 0 ; i_FMM < Nd ;  i_FMM++){ 
	V[i_FMM].Type = SCALAR ;
	V[i_FMM].Val[0] = V[i_FMM].Val[MAX_DIM] = 0. ;
	kr = (xxs*Current.FMM.Kdir[i_FMM][0] + 
	      yys*Current.FMM.Kdir[i_FMM][1] +
	      zzs*Current.FMM.Kdir[i_FMM][2])/r ; 

	LegendreRecursive( ns, 0, kr, P);
    	
	for (j_FMM = 0 ; j_FMM <= ns ; j_FMM++){ 
	  s = ((j_FMM/2)%2)? -1 : 1 ;	 
	  sbjn = jsph[j_FMM] ;
	  sbyn = ysph[j_FMM] ;
	  Pl = (2*j_FMM+1)* P[j_FMM] ;

	  V[i_FMM].Val[0]       += s * Pl * ((j_FMM%2) ?   sbjn : sbyn ) ;
	  V[i_FMM].Val[MAX_DIM] += s * Pl * ((j_FMM%2) ? (-sbyn): sbjn ) ;
	}

	V[i_FMM].Val[0]       *= cr*Current.FMM.Weight[i_FMM];
	V[i_FMM].Val[MAX_DIM] *= cr*Current.FMM.Weight[i_FMM];
      }

      break;      
    default :
      Msg(GERROR, "Case 3D: Bad Flag_GF 'GF_HelmholtzThinWire' (%d)", Current.FMM.Flag_GF);  
      break;

    }
    break;

  default :
    Msg(GERROR, "Bad Parameter for 'GF_HelmholtzThinWire' (%d)", (int)Fct->Para[0]);
    break;
  }
     
  GetDP_End ;
}




/* ------------------------------------------------------------------------ */
/*  G F _ G r a d H e l m h o l t z                                         */
/* ------------------------------------------------------------------------ */

/* the gradient is taken relative to the destination point (x,y,z) */

void GF_GradHelmholtz (F_ARG) {