gf_laplacexform.c

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

#define RCSID "$Id: GF_LaplacexForm.c,v 1.19 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 "BF_Function.h"
#include "CurrentData.h"
#include "Data_DefineE.h"
#include "Numeric.h"
#include "Numeric_F.h"


#define F_ARG2						\
  struct Element * Element, struct Function * Fct,	\
  void  (*xFunctionBF) (), int NumEntity, 		\
  double x, double y, double z, struct Value * Val

#define CAST  void(*)()

#define MAX_NODES        6 
#define EPSILON          1.e-8
#define EPSILON2         1.e-20

/* this is a hack... */
#define RADIUS           0.154797 


/* ------------------------------------------------------------------------ */
/*  G F _ L a p l a c e x F o r m                                           */
/* ------------------------------------------------------------------------ */

void GF_LaplacexForm (F_ARG2) {
  
  double   xs[MAX_NODES], ys[MAX_NODES], zs[MAX_NODES], u[3], v[3], n[3];
  double   u2=0., v2=0., xl=0., yl=0., zl=0., zl_2=0. ;
  double   Area, m0[3], m1[3], m2[3] ;
  int      Type_Int=0, i, j = 1 ;
  double   a=0., b=0., c=0., d, e, f, i1, I1 = 0., Iua, Iva, r2;
  double   s0m=0., s0p=0., s1m=0., s1p=0., s2m=0., s2p=0., t00, t10, t20, t0m_2, t0p_2, t1p_2;
  double   r00_2=0., r10_2=0., r20_2=0., r00, r10, r20, r0p=0., r0m=0., r1p=0.;
  double   f20=0., f21=0., f22=0., B0, B1, B2 ;
  double   f30, f31, f32, N10, N20, N30 ;
  double   DetJ, valr, vali ;

  GetDP_Begin("GF_LaplacexForm");

  Val->Val[MAX_DIM] = 0.0 ;
  
  switch ((int)Fct->Para[0]) {
    
  case _2D : 
    
    switch (Element->ElementSource->Type) {
      
    case POINT :      
      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;

      r2 = SQU(x-xs[0])+SQU(y-ys[0]) ;
      if (r2 > SQU(RADIUS)){
	Val->Type = SCALAR ;
	Val->Val[0] = - ONE_OVER_FOUR_PI * log(r2) ;
      }
      else{
	Val->Type = SCALAR ;
	Val->Val[0] = - ONE_OVER_FOUR_PI * log(SQU(RADIUS)) ;
      }
      break ;
     
    case LINE :
      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      
      if(xFunctionBF == (CAST)BF_Volume) {      
	a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) ;
	b = 2. * ((x-xs[0])*(xs[0]-xs[1]) + (y-ys[0])*(ys[0]-ys[1])) ;
	c = SQU(x-xs[0]) + SQU(y-ys[0]) ;     
	d = 0.5 * b / a ;  
	e = c / a ;  
	f = e - d*d ;
	
	if      (f > EPSILON)      { Type_Int = 1; }
	else if (fabs(f) < EPSILON){ Type_Int = 0; }
	else                       { Type_Int = -1; f = -f; }	
	if (Element->Num == Element->ElementSource->Num)  Type_Int = 2 ;
	if ((c == 0) || ((b == -2*a) && (c == a)))  Type_Int = 3 ;
	
	switch (Type_Int) {   
	case -1 :	
	  I1 = log(a) + 
	    ( (d+1.) * log(SQU(d+1.) - f) - 2.*(d+1.) + 
	      sqrt(f) * log((d+1.+sqrt(f))/(d+1.-sqrt(f))) ) -
	    ( d*log(d*d-f) - 2.*d + sqrt(f)*log((d+sqrt(f))/(d-sqrt(f))) ) ;
	  break ;
	case 0 :
	  I1 = log(a) + (d+1.)*log(SQU(d+1.)) - d*log(SQU(d)) - 2. ; 
	  break ;      	    
	case 1 :
	  I1 = log(a) + 
	    ( (d+1.) * log(SQU(d+1.) + f) - 2.*(d+1.) + 
	      2.*sqrt(f) * atan((d+1.)/sqrt(f)) ) -
	    ( d*log(d*d+f) - 2.*d + 2.*sqrt(f)*atan(d/sqrt(f)) )  ;
	  break ;	  	    
	case 2 : 
	  i1 = -b / (2.*a) ;
	  I1 = 2. * i1 * (log(i1) - 1.) + 
	    2. * (1.-i1) * (log(1.-i1) - 1.) + log(a) ;
	  break ;	    
	case 3 : 
	  I1 = .5 * log(a) - 1. ;
	  break ;    
	} 
	
	Val->Type = SCALAR ;
	Val->Val[0] = - ONE_OVER_FOUR_PI * I1 ; 
      }
      else {
	Msg(GERROR, "Unknown Basis Function Type for 'GF_LaplacexForm'");
      }
      break ;
      
      
    case TRIANGLE :
      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;

      if(xFunctionBF == NULL) {
	clt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);	
	DetJ = (xs[2]-xs[0]) * (ys[1]-ys[0]) - (xs[1]-xs[0]) * (ys[2]-ys[0]) ;
	Val->Type = SCALAR ;
	Val->Val[0] = valr * DetJ ;
      }
      else if(xFunctionBF == (CAST)BF_Volume) {
	clt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	Val->Type = SCALAR ;
	Val->Val[0] = valr * 2  /*  *DetJ/DetJ  */ ; 
      }
      else if(xFunctionBF == (CAST)BF_Node) {
	switch(NumEntity){
	case 1 : clt2dl_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali); break;
	case 2 : clt2dl_(&x,&y,&xs[1],&xs[2],&xs[0],&ys[1],&ys[2],&ys[0],&valr,&vali); break;
	case 3 : clt2dl_(&x,&y,&xs[2],&xs[0],&xs[1],&ys[2],&ys[0],&ys[1],&valr,&vali); break;
	}
	DetJ = (xs[2]-xs[0]) * (ys[1]-ys[0]) - (xs[1]-xs[0]) * (ys[2]-ys[0]) ;
	Val->Type = SCALAR ;
	Val->Val[0] = valr * DetJ ; 
      }
      else{
	Msg(GERROR, "Unknown Basis Function Type for 'GF_LaplacexForm'");
      }
      break ;


    case QUADRANGLE :
      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;
      xs[3] = Element->ElementSource->x[3] ; ys[3] = Element->ElementSource->y[3] ;
      
      if(xFunctionBF == NULL) {
	clt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	DetJ = (xs[2]-xs[0]) * (ys[1]-ys[0]) - (xs[1]-xs[0]) * (ys[2]-ys[0]) ;
	Val->Val[0] = valr * DetJ ; 
	clt2d_(&x,&y,&xs[0],&xs[2],&xs[3],&ys[0],&ys[2],&ys[3],&valr,&vali);
	DetJ = (xs[3]-xs[0]) * (ys[2]-ys[0]) - (xs[2]-xs[0]) * (ys[3]-ys[0]) ;
	Val->Val[0] += valr * DetJ ;       
	Val->Type = SCALAR ;
      }
      else if(xFunctionBF == (CAST)BF_Volume) {
	clt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	Val->Val[0] = valr * 2  /*  *DetJ/DetJ  */ ; 
	clt2d_(&x,&y,&xs[0],&xs[2],&xs[3],&ys[0],&ys[2],&ys[3],&valr,&vali);
	Val->Val[0] += valr * 2  /*  *DetJ/DetJ  */ ; 
	Val->Type = SCALAR ;
      }
      else{
	Msg(GERROR, "Unknown Basis Function Type for 'GF_LaplacexForm'");
      }
      break ;
      
    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_LaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));
    }
    
    break;
    


  case _3D :
    
    switch (Element->ElementSource->Type) {
      
    case LINE :
      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      zs[0] = Element->ElementSource->z[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      zs[1] = Element->ElementSource->z[1] ;
      
      a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) + SQU(zs[0]-zs[1]) ;
      b = 2. * ((x-xs[0])*(xs[0]-xs[1]) + 
		(y-ys[0])*(ys[0]-ys[1]) + 
		(z-zs[0])*(zs[0]-zs[1])) ;
      c = SQU(x-xs[0]) + SQU(y-ys[0]) + SQU(z-zs[0]) + SQU(RADIUS) ;
      
      Val->Val[0] = ONE_OVER_FOUR_PI *
	log( ( 2.*sqrt(a*(a+b+c))+2.*a+b ) / ( 2.*sqrt(a*c)+b ) ) ;
      Val->Type = SCALAR ;
      break ;
    
    case TRIANGLE :
    case QUADRANGLE :
      if(xFunctionBF == (CAST)BF_Volume) Type_Int = 1 ;
      if(xFunctionBF == (CAST)BF_Node)   Type_Int = 2 ;

      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      zs[0] = Element->ElementSource->z[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      zs[1] = Element->ElementSource->z[1] ;
      xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;
      zs[2] = Element->ElementSource->z[2] ;

      if (Element->ElementSource->Type ==  QUADRANGLE) {
	     xs[3] = Element->ElementSource->x[3] ; ys[3] = Element->ElementSource->y[3] ;
	     zs[3] = Element->ElementSource->z[3] ;
	     j = 0 ;
       };
      
      for(i = j; i < 2; i++){      
	/* triangle side lengths */
	a = sqrt(SQU(xs[1]-xs[0]) + SQU(ys[1]-ys[0]) +  SQU(zs[1]-zs[0]));
	b = sqrt(SQU(xs[2]-xs[1]) + SQU(ys[2]-ys[1]) +  SQU(zs[2]-zs[1]));
	c = sqrt(SQU(xs[2]-xs[0]) + SQU(ys[2]-ys[0]) +  SQU(zs[2]-zs[0]));

	/* local system (u,v,w) centered at (xs[0],ys[0],zs[0]) */
	u[0] = (xs[1]-xs[0])/a;
	u[1] = (ys[1]-ys[0])/a;
	u[2] = (zs[1]-zs[0])/a;
      
	/* triangle normal */
	Geo_CreateNormal(Element->ElementSource->Type,xs,ys,zs,n);
	
	/* v = n /\ u */ 
	v[0] = n[1]*u[2]-n[2]*u[1];
	v[1] = n[2]*u[0]-n[0]*u[2];
	v[2] = n[0]*u[1]-n[1]*u[0];
	
	u2 = (xs[2]-xs[0])*u[0] + (ys[2]-ys[0])*u[1] + (zs[2]-zs[0])*u[2]; /* u2 coordinate */
	v2 = (xs[2]-xs[0])*v[0] + (ys[2]-ys[0])*v[1] + (zs[2]-zs[0])*v[2];  /* triangle height, v2 coordinate */
     
	/* local coordinates of the observation point (xl, yl, zl) */ 
	xl = u[0] * (x-xs[0]) +  u[1] * (y-ys[0]) + u[2] * (z-zs[0]);
	yl = v[0] * (x-xs[0]) +  v[1] * (y-ys[0]) + v[2] * (z-zs[0]);
	zl = n[0] * (x-xs[0]) +  n[1] * (y-ys[0]) + n[2] * (z-zs[0]);
		
	s0m = -( (a-xl) * (a-u2) + yl*v2 ) / b;
	s0p = s0m + b;
	s1p = ( xl * u2 + yl * v2 ) / c;
	s1m = s1p - c;
	s2m = - xl;
	s2p = a - xl;
	
	/* distance observation point projection on triangle plane to triangle local vertices*/
	/* t1m = t0p ; t2p = t0m ; t2m = t1p ; */ 
	t00 = (yl * (u2-a) + v2 * (a-xl)) / b;
	t10 = (xl * v2 - yl * u2) / c;
	t20 = yl;

	t0m_2 = (a-xl)*(a-xl) + yl*yl; 
	t0p_2 = (u2-xl)*(u2-xl) + (v2-yl)*(v2-yl); 
	t1p_2 = xl*xl + yl*yl;
	
	/* minimum distances^2 from the observation point to each triangle side*/
	zl_2 = SQU(zl) ;
	r00_2 = SQU(t00) + zl_2 ;
	r10_2 = SQU(t10) + zl_2 ;
	r20_2 = SQU(t20) + zl_2 ;
	
	/* distances from observation point to the vertices*/
	r0p = sqrt(t0p_2 + zl_2);      
	r0m = sqrt(t0m_2 + zl_2);
	r1p = sqrt(t1p_2 + zl_2);

	r00 = sqrt(r00_2);
	r10 = sqrt(r10_2);
	r20 = sqrt(r20_2);
	
	/* intermediate functions */ 
	if(r00 <= EPSILON*(fabs(s0m)+fabs(s0p)) ){ f20 = log(s0m/s0p) ; B0 = 0; }
	else{
	  if (!(r0m + s0m)) 
	    Msg(GERROR,"1/0 in GF_LaplacexForm (case _3D TRIANGLE) Num %d Obs %.15e %.15e %.15e",
		Element->ElementSource->Num, x, y, z) ;
	  f20 = log((r0p + s0p) / (r0m + s0m));
	  B0  = atan(t00*s0p/(r00_2+fabs(zl)*r0p))-atan(t00*s0m/(r00_2+fabs(zl)*r0m));
	}
	
	if(r10 <= EPSILON*(fabs(s1m)+fabs(s1p)) ){ f21 = log(s1m/s1p); B1 = 0; }
	else{
	  if(!(r0p + s1m)) 
	    Msg(GERROR,"1/0 in GF_LaplacexForm (case _3D TRIANGLE) Num %d Obs %.15e %.15e %.15e",
		Element->ElementSource->Num, x, y, z) ;
	  f21 = log((r1p + s1p) / (r0p + s1m));
	  B1 =  atan(t10*s1p/(r10_2+fabs(zl)*r1p))-atan(t10*s1m/(r10_2+fabs(zl)*r0p));
	}
	
	if(r20 <= EPSILON*(fabs(s2m)+fabs(s2p)) ){ f22 = log(s2m/s2p); B2 = 0; }
	else{
	  if(!(r1p+s2m))
	    Msg(GERROR,"1/0 in GF_LaplacexForm (case _3D TRIANGLE) Num %d Obs %.15e %.15e %.15e",
		Element->ElementSource->Num, x, y, z) ;
	  f22 = log((r0m + s2p) / (r1p + s2m));
	  B2 = atan(t20*s2p/(r20_2+fabs(zl)*r0m))-atan(t20*s2m/(r20_2+fabs(zl)*r1p));
	}
			
	I1 += -fabs(zl)*(B0+B1+B2) + t00*f20+t10*f21+t20*f22 ; /* 1/r integral solution*/
	
	if (j == 0){ xs[1] = xs[2]; ys[1] = ys[2]; zs[1] = zs[2]; xs[2] = xs[3]; ys[2] = ys[3]; zs[2] = zs[3];}
      }

      
      switch ( Type_Int ){
      case 1 : /* BF_Volume */
	Area = a * v2/2 ;/* Triangle area */
	Val->Val[0] = I1 /Area ; 	
	break;

      case 2 : /* BF_Node */
	if (!v2) 	
	  Msg(GERROR,"1/0 in GF_LaplacexForm (case _3D TRIANGLE) v2 %e", v2);

	f30 = (s0p*r0p-s0m*r0m) + r00_2 * f20 ; /* f3i */
	f31 = (s1p*r1p-s1m*r0p) + r10_2 * f21 ;
	f32 = (s2p*r0m-s2m*r1p) + r20_2 * f22 ;

	m0[0] = ((ys[2] - ys[1]) * n[2] - (zs[2] - zs[1]) * n[1])*f30/b ;
	m0[1] = ((zs[2] - zs[1]) * n[0] - (xs[2] - xs[1]) * n[2])*f30/b ;
	m0[2] = ((xs[2] - xs[1]) * n[1] - (ys[2] - ys[1]) * n[0])*f30/b ;
	
	m1[0] = ((ys[0] - ys[2]) * n[2] - (zs[0] - zs[2]) * n[1])*f31/c ;
	m1[1] = ((zs[0] - zs[2]) * n[0] - (xs[0] - xs[2]) * n[2])*f31/c ;
	m1[2] = ((xs[0] - xs[2]) * n[1] - (ys[0] - ys[2]) * n[0])*f31/c ;
	
	m2[0] = (u[1] * n[2] - u[2]* n[1])*f32 ;
	m2[1] = (u[2] * n[0] - u[0]* n[2])*f32 ;
	m2[2] = (u[0] * n[1] - u[1]* n[0])*f32 ;
	
	
	Iua = (u[0] * (m0[0] + m1[0] + m2[0]) +
               u[1] * (m0[1] + m1[1] + m2[1]) + 
	       u[2] * (m0[2] + m1[2] + m2[2]))/2 ; 
	
	Iva = (v[0] * (m0[0] + m1[0] + m2[0]) +
               v[1] * (m0[1] + m1[1] + m2[1]) + 
	       v[2] * (m0[2] + m1[2] + m2[2]))/2 ; 
	
	switch(NumEntity){
	case 1 :  
	  N10 = 1 - xl/a + (u2/a -1) * yl/v2 ;
	  Val->Val[0] = N10 * I1  - Iua/a + (u2/a-1) * Iva/v2 ;
	  break;
	case 2 : 
 	  N20 = xl/a - u2/a * yl/v2 ;
 	  Val->Val[0] = N20 * I1 + Iua/a - u2/a * Iva/v2 ;
	  break;
	case 3 :  
	  N30 = yl/v2 ;
	  Val->Val[0] = N30 * I1 + Iva/v2 ;
	  break;
	}
	break;
      default :
	Msg(GERROR, "Unknown Basis Function Type for 'GF_LaplacexForm'");
      }

      Val->Val[0] *= ONE_OVER_FOUR_PI ;
      if (j == 0){ Val->Val[0] /= 2; }
      Val->Type = SCALAR ;
      break ;
    
    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_LaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));
    }
    break ;

  default :
    Msg(GERROR, "Unknown Dimension (%d) for 'GF_LaplacexForm'", 
	(int)Fct->Para[0]);
    
  }
  
  GetDP_End ;
}



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

void GF_GradLaplacexForm (F_ARG2) {
  
  double  xs[MAX_NODES], ys[MAX_NODES], zs[MAX_NODES] ;
  double  xxs, yys, r2, EPS ;
  double  a, b, c, a2, I1, I2 ;

  double  mx=0., my=0., valr, vali, DetJ ;
  double  f0[3], f1[3], f2[3], N10, N20, N30 ;
  double  m0[3], m1[3], m2[3], s0[3], s1[3] ;
  double  umf2i, us0, us1, us2, vmf2i, vs0, vs1, vs2 ;
  double  u[3], v[3], n[3], u2, v2, xl, yl, zl, zl_2 ;
  double  area, I[3], Iua[3], Iva[3] ;
  double  s0m, s0p, s1m, s1p, s2m, s2p, t00, t10, t20, t0m_2, t0p_2, t1p_2;
  double  r00_2, r10_2, r20_2,  r00, r10, r20, r0p, r0m, r1p, f20, f21, f22, B0, B1, B2, B ;
  int Type_Int=0 ;

  GetDP_Begin("GF_GradLaplacexForm");

  Val->Val[MAX_DIM] =  Val->Val[MAX_DIM + 1] = Val->Val[MAX_DIM + 2] = 0. ;

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