gf_laplace.c

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

#define RCSID "$Id: GF_Laplace.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 "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 
*/

/* ------------------------------------------------------------------------ */
/*  G F _ L a p l a c e                                                     */
/* ------------------------------------------------------------------------ */

extern int Flag_RemoveSingularity ;

void  GF_Laplace (F_ARG) {

  double d, cte, r, phi_r, r_l, r_l_1, F, a, b, c1, c2  ;
  double xxs, yys, zzs, Theta, Phi, cTheta ;
  double Plm ; 

  int i, k, l, m, Nd ; 

  GetDP_Begin("GF_Laplace");


  switch((int)Fct->Para[0]){

  case _1D : /* r/2 */
    V->Val[0] = 0.5 * sqrt(SQU(Current.x-Current.xs)) ;  
    break;

  case _2D : /* 1/(2*Pi) * ln(1/r) = -1/(4*Pi)*ln(r^2) */
    switch(Current.FMM.Flag_GF){
    case FMM_DIRECT :
      d = SQU(Current.x-Current.xs) + SQU(Current.y-Current.ys) ;
      if(!d) Msg(GERROR, "Log(0) in 'GF_Laplace'") ;
      V->Val[0] = - ONE_OVER_FOUR_PI * log(d) ;
      V->Val[MAX_DIM] = 0. ; 
      break;

    case FMM_AGGREGATION :
      Nd = Current.FMM.NbrDir ;
      /* Normalization with regard to the maximum distance in the FMM group, (x-xc)/Rx */
      r = sqrt(SQU(Current.xs-Current.FMM.Xgc) + SQU(Current.ys-Current.FMM.Ygc))/Current.FMM.Rsrc ; 
      phi_r = atan2(Current.ys-Current.FMM.Ygc, Current.xs-Current.FMM.Xgc) ;    
      
      cte = - ONE_OVER_TWO_PI ;
      V[0].Type = SCALAR ;
      V[0].Val[0]       = cte ;
      V[0].Val[MAX_DIM] = 0. ;
      
      for (i = 1 ; i < Nd ; i++){
	cte *= r ; 
	V[i].Type = SCALAR ;
	V[i].Val[0]       = cte * cos(i*phi_r) ;
	V[i].Val[MAX_DIM] = cte * sin(i*phi_r) ;
      }
      break;

    case FMM_DISAGGREGATION :
     Nd = Current.FMM.NbrDir ;
     /* Normalization with regard to the maximum distance in the FMM group,   (yc-y)/Ry */
     r = sqrt(SQU(Current.FMM.Xgc-Current.x) + SQU(Current.FMM.Ygc-Current.y))/Current.FMM.Robs ;
     phi_r = atan2(Current.FMM.Ygc-Current.y, Current.FMM.Xgc-Current.x) ;
     cte = 1. ;
     V[0].Type = SCALAR ;
     V[0].Val[0]       = cte ;
     V[0].Val[MAX_DIM] = 0.  ;
     for (i = 1 ; i < Nd ;  i++){ 
       cte *= r ;
       V[i].Type = SCALAR ;
       V[i].Val[0]       = cte * cos(i*phi_r) ;
       V[i].Val[MAX_DIM] = cte * sin(i*phi_r) ;
     }
     break ;

    case FMM_TRANSLATION : /* 2D */
      Nd = Current.FMM.NbrDir ;/* (x, y, z) Destination, (xs, ys, zs) Origin */
      r = sqrt( SQU(Current.x-Current.xs) + SQU(Current.y-Current.ys)) ; /* (yc-xc) */  
      phi_r = atan2(Current.y-Current.ys, Current.x-Current.xs) ;

      if(!r) Msg(GERROR, "1/0 in 'GF_Laplace (Translation - 2D)'") ;
     
      a = 1. ;  c1 = 1. ;
      for (i = 0 ; i < Nd ;  i++){
	b = 1. ; c2 = 1. ;
	for (k = 0 ; k < Nd ; k++){
	  if ( i == 0  && k == 0){
	    V[0].Type = SCALAR ;
	    V[0].Val[0]       = log(r) ; /* log(yc-xc) */
	    V[0].Val[MAX_DIM] = phi_r  ;
	  }
	  else{
	    F = -BinomialCoef(i+k, k)/(i+k);
	    /* Normalization (Rx^i * Ry^k) */
	    cte  = F * a * b * c1 *c2 ;
	    V[i*Nd + k].Type = SCALAR ;
	    V[i*Nd + k].Val[0]       =   cte * cos((i+k)*phi_r) ;
	    V[i*Nd + k].Val[MAX_DIM] = - cte * sin((i+k)*phi_r) ;
	  }      
	  b *= Current.FMM.Rsrc;
	  c2 /= r;
	}
	a *= Current.FMM.Robs;
	c1 /= r ;
      }       
      break;
 
   default :
      Msg(GERROR, "Case 2D: Bad Flag_GF 'GF_Laplace' (%d)", Current.FMM.Flag_GF);  
      break;
    }
    break;

  case _3D : /* 1/(4*Pi*r) */
    switch(Current.FMM.Flag_GF){          
    case FMM_DIRECT :         
      if(0){/* Flag_RemoveSingularity */
	printf("  ->removed in %d %d \n", Current.Element->Num, Current.ElementSource->Num);
	V->Val[0] = ONE_OVER_FOUR_PI ;
      }
      else{
	d = SQU(Current.x-Current.xs) + SQU(Current.y-Current.ys) + SQU(Current.z-Current.zs) ;
	if(!d) Msg(GERROR, "1/0 in 'GF_Laplace'") ;
        V->Val[0] = ONE_OVER_FOUR_PI / sqrt(d) ;
      }
      V->Val[MAX_DIM] = 0. ; 
      break;
      /* Not normalization for 3D */
    case FMM_AGGREGATION : /* 3D */
      Nd = Current.FMM.NbrDir ;
      xxs = Current.xs - Current.FMM.Xgc ;
      yys = Current.ys - Current.FMM.Ygc ;
      zzs = Current.zs - Current.FMM.Zgc ;

      r = sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)) ;
      Theta = atan2(sqrt(SQU(xxs)+SQU(yys)), zzs) ; 
      Phi   = atan2(yys, xxs) ;
      cTheta = cos(Theta);

      r_l = 1.;
      for (l = 0 ; l < Nd ; l++){
	for (m = -l ; m <= l ; m++){
	  Plm  = Legendre(l, m, cTheta) ;	 	  
	  cte = r_l/Factorial(l+m) * Plm ;  
	  V[l*(l+1)+m].Type = SCALAR ;  	  
	  V[l*(l+1)+m].Val[0]       =  cte * cos(m*Phi) ;
	  V[l*(l+1)+m].Val[MAX_DIM] = -cte * sin(m*Phi) ;
	}
	r_l *= r ;
      }
      break;

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

      r = sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)) ;  
      Theta = atan2(sqrt(SQU(xxs)+SQU(yys)), zzs) ;
      Phi   = atan2(yys, xxs) ;      
      cTheta = cos(Theta) ;

      r_l = 1. ;
      for (l = 0 ; l < Nd ; l++){
	for (m = -l ; m <= l ; m++){
	  Plm  = Legendre(l, m, cTheta) ;
	  cte = r_l/Factorial(l+m) * Plm ;

	  V[l*(l+1)+m].Type = SCALAR ;  	  
	  V[l*(l+1)+m].Val[0]       =  cte * cos(m*Phi) ;
	  V[l*(l+1)+m].Val[MAX_DIM] = -cte * sin(m*Phi) ;
	}
	r_l *= r ;
      }
      break;

    case FMM_TRANSLATION: /* 3D */
      Nd =  Current.FMM.NbrDir ;     
      xxs = Current.x - Current.xs ;
      yys = Current.y - Current.ys ;
      zzs = Current.z - Current.zs ;
      
      r =  sqrt(SQU(xxs) + SQU(yys) + SQU(zzs)); /* Distance between FMM group centers */ 
      if(!r) Msg(GERROR, "1/r in 'GF_Laplace (Translation - 3D)'") ;

      Theta = atan2(sqrt(SQU(xxs)+SQU(yys)), zzs) ; 
      Phi   = atan2(yys, xxs) ;      
      cTheta = cos(Theta) ;

      r_l_1 = ONE_OVER_FOUR_PI ;
      for (l = 0 ; l < 2*Nd ; l++){
	r_l_1 /= r ;
	for (m = -l ; m <= l ; m++){
	  Plm  = Legendre(l, m, cTheta) ;  	  
	  cte = Factorial((double)(l-m)) * r_l_1 * Plm ;
	  V[l*(l+1)+m].Type = SCALAR ; 
	  V[l*(l+1)+m].Val[0]       =  cte * cos(m*Phi) ;
	  V[l*(l+1)+m].Val[MAX_DIM] =  cte * sin(m*Phi) ;
	}
      }	
      break;

    default :
      Msg(GERROR, "Case 3D: Bad Flag_GF 'GF_Laplace' (%d)", Current.FMM.Flag_GF);  
      break;
    }    
    break;
    
  default :
    Msg(GERROR, "Bad Parameter for 'GF_Laplace' (%d)", (int)Fct->Para[0]);
    break;
  }

  V->Type = SCALAR ;
  V->Val[MAX_DIM] = 0. ;  

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G F _ G r a d L a p l a c e                                             */
/* ------------------------------------------------------------------------ */

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

void GF_GradLaplace (F_ARG) {
  double xxs, yys, zzs, Theta, Phi, sTheta, cTheta, sPhi, cPhi, Plm, dPlm;
  double cte, r, phi_r, cphi_r, sphi_r, cr, r_l_2, a, b ;
  double r_l, r_l_1, r_2, rxy, rxy_2, srxy, c1, tx1, ty1, tx2, ty2 ;
  int i, Nd, k, l, m ;   

  GetDP_Begin("GF_GradLaplace");

  V->Type = VECTOR ;

  switch((int)Fct->Para[0]){
  case _2D :
    switch(Current.FMM.Flag_GF){
    case FMM_DIRECT :
      xxs = Current.x-Current.xs ; 
      yys = Current.y-Current.ys ; 
      r = SQU(xxs)+SQU(yys) ;
      if(!r) Msg(GERROR, "1/0 in 'GF_GradLaplace'") ;
     
      V->Val[0] = - ONE_OVER_TWO_PI * xxs / r ;
      V->Val[1] = - ONE_OVER_TWO_PI * yys / r ;
      V->Val[2] = 0.0 ;
      V->Val[MAX_DIM    ] = V->Val[MAX_DIM + 1] = V->Val[MAX_DIM + 2] = 0. ;
      break ;    


    case FMM_AGGREGATION : /* 2D -> SCALAR */
      Nd = Current.FMM.NbrDir ;
      r = sqrt(SQU(Current.FMM.Xgc-Current.xs) + SQU(Current.FMM.Ygc-Current.ys))/Current.FMM.Rsrc;  /* (x-xc)/Rx */
      phi_r = atan2( Current.FMM.Ygc-Current.ys, Current.FMM.Xgc-Current.xs) ;    
      cte =  - ONE_OVER_TWO_PI/r ; 
      for (i = 0 ; i < Nd ; i++){
	cte *= r; 
	V[i].Type = SCALAR ;
	V[i].Val[0]       = cte * cos(i*phi_r) ;
	V[i].Val[MAX_DIM] = cte * sin(i*phi_r) ;
      } 
      break;

    case FMM_DISAGGREGATION : /* 2D -> VECTOR */
     Nd = Current.FMM.NbrDir ;
     xxs = -Current.FMM.Xgc + Current.x;
     yys = -Current.FMM.Ygc + Current.y;
     r = sqrt(SQU(xxs) + SQU(yys))/Current.FMM.Robs;
     phi_r = atan2(yys, xxs) ;
     if(!r) Msg(GERROR, "1/0 in 'GF_GradLaplace (Disaggregation - 2D)'") ;
     cte = -1/r/r ; 
     for (i = 0 ; i < Nd ; i++){
	cte *= r ; 
	cphi_r = cos((i+1)*phi_r) ;
	sphi_r = sin((i+1)*phi_r) ;
	V[i].Type       = VECTOR ;
	V[i].Val[0]         = (xxs * cphi_r + yys * sphi_r) * cte;
	V[i].Val[1]         = (yys * cphi_r - xxs * sphi_r) * cte;
	V[i].Val[2]         = 0.;
	V[i].Val[MAX_DIM]   = -V[i].Val[1];	
	V[i].Val[MAX_DIM+1] =  V[i].Val[0];
	V[i].Val[MAX_DIM+2] = 0.;
      }
     break ;

    case FMM_TRANSLATION : /* 2D */
      Nd = Current.FMM.NbrDir ;  /* (x, y, z) Destination, (xs, ys, zs) Origin */ 
      r = sqrt(SQU(Current.xs-Current.x) + SQU(Current.ys-Current.y)) ;  
      phi_r = atan2(Current.ys-Current.y, Current.xs-Current.x) ;
      
      if(!r) Msg(GERROR, "1/0 in 'GF_GradLaplace (Translation - 2D)'") ;

      a = 1/Current.FMM.Rsrc ;
      for (i = 0 ; i < Nd ;  i++){
	a *= Current.FMM.Rsrc/r ;
	b = r/Current.FMM.Robs/Current.FMM.Robs ;	
	for (k = 0 ; k < Nd ; k++){
	  b *= Current.FMM.Robs/r ;
	  cte = BinomialCoef(i+k, k)* a * b  ;

	  V[i*Nd + k].Type = SCALAR ;
	  V[i*Nd + k].Val[0]       =   cte * cos((i+k+1)*phi_r) ;
	  V[i*Nd + k].Val[MAX_DIM] = - cte * sin((i+k+1)*phi_r) ;
 	}	
      }        

      break;

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

  case _3D :
    switch(Current.FMM.Flag_GF){
    case FMM_DIRECT :
      xxs = Current.x-Current.xs ; 
      yys = Current.y-Current.ys ; 
      zzs = Current.z-Current.zs ;
      r = CUB(sqrt(SQU(xxs)+SQU(yys)+SQU(zzs))) ;
      
      if(!r) Msg(GERROR, "1/0 in 'GF_GradLaplace'") ;
     
      V->Val[0] = - ONE_OVER_FOUR_PI * xxs / r ;
      V->Val[1] = - ONE_OVER_FOUR_PI * yys / r ;
      V->Val[2] = - ONE_OVER_FOUR_PI * zzs / r ;
      V->Val[MAX_DIM    ] = V->Val[MAX_DIM + 1] = V->Val[MAX_DIM + 2] = 0. ;
      break ;

    case FMM_AGGREGATION : /* 3D */
      Nd = Current.FMM.NbrDir ;

      xxs = Current.xs - Current.FMM.Xgc ;
      yys = Current.ys - Current.FMM.Ygc ;
      zzs = Current.zs - Current.FMM.Zgc ;

      r = sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)) ;
      Theta = atan2( sqrt(SQU(xxs)+SQU(yys)), zzs ) ; 
      Phi   = atan2(yys, xxs) ;
      cTheta = cos(Theta);

      r_l = 1. ;
      for (l = 0 ; l < Nd ; l++){
	for (m = -l ; m <= l ; m++){
	  Plm  = Legendre(l, m, cTheta) ;
	  cte = r_l/Factorial(l+m) * Plm ;

	  V[l*(l+1)+m].Type = SCALAR ;  	  
	  V[l*(l+1)+m].Val[0]       =  cte * cos(m*Phi) ;
	  V[l*(l+1)+m].Val[MAX_DIM] = -cte * sin(m*Phi) ;
	}
	r_l *= r ; 
      }
      break;

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

      r_2   = SQU(xxs)+SQU(yys)+SQU(zzs) ;
      r = sqrt(r_2) ;
      if(!r) Msg(GERROR, "1/0 in 'GF_GradLaplace (Disaggregation - 3D)'") ;

      rxy_2 = SQU(xxs)+SQU(yys) ;
      rxy =  sqrt(rxy_2) ;

      Theta = atan2(rxy, zzs) ; 
      Phi   = atan2(yys, xxs) ;
      cTheta = cos(Theta) ;
      sTheta = sin(Theta) ;
      srxy =  sTheta * rxy ;

      c1 = sTheta *  zzs/rxy ;