Skip to content
Snippets Groups Projects
Select Git revision
  • 03c4ceb15fe301b9c0fae88aa8f2ad48e9d34dff
  • master default protected
  • dof-renumbering
  • test-dof-hash
  • gdemesy-master-patch-30528
  • eval-space-time
  • oscillating_multiharm
  • MH_movement
  • axisqu
  • write_vtu_and_ensight_formats
  • movingband
  • CP_1972_add_vtu_file_writing
  • mortar
  • fast_freq_sweep_Resolution
  • applyresolvent_again
  • marteaua-master-patch-54323
  • patch-1
  • binde-master-patch-08072
  • binde-master-patch-52461
  • BCGSL
  • resolvent
  • getdp_3_5_0
  • getdp_3_4_0
  • getdp_3_3_0
  • getdp_3_2_0
  • getdp_3_1_0
  • getdp_3_0_4
  • getdp_3_0_3
  • getdp_3_0_2
  • getdp_3_0_1
  • getdp_3_0_0
  • onelab_mobile_2.1.0
  • getdp_2_11_3 protected
  • getdp_2_11_2 protected
  • getdp_2_11_1 protected
  • getdp_2_11_0 protected
  • getdp_2_10_0 protected
  • getdp_2_9_2 protected
  • getdp_2_9_1 protected
  • getdp_2_9_0 protected
  • getdp_2_8_0 protected
41 results

GF_HelmholtzxForm.cpp

Blame
  • GF_HelmholtzxForm.cpp 5.64 KiB
    // GetDP - Copyright (C) 1997-2018 P. Dular and C. Geuzaine, University of Liege
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <getdp@onelab.info>.
    //
    // Contributor(s):
    //   Ruth Sabariego
    //
    
    #include <math.h>
    #include "ProData.h"
    #include "ProDefine.h"
    #include "BF.h"
    #include "GF.h"
    #include "Gauss.h"
    #include "GeoData.h"
    #include "Message.h"
    
    #define SQU(a)     ((a)*(a)) 
    #define CUB(a)     ((a)*(a)*(a)) 
    #define ONE_OVER_TWO_PI    1.5915494309189534E-01
    
    #define MAX_NODES  6 
    #define EPSILON    1.e-10
    
    /* ------------------------------------------------------------------------ */
    /*  G F _ H e l m h o l t z x F o r m                                       */
    /* ------------------------------------------------------------------------ */
    
    void GF_HelmholtzxForm(GF_ARGX)
    {
      double   xs[MAX_NODES], ys[MAX_NODES], zs[MAX_NODES], u[3], v[3], n[3], l2[2];
      double   xl, yl, zl ;
      int      i, j = 1 ;
      int      i_IntPoint, NGT_Points = 7;
      double   a, b, c, d, us, vs, ws, usi, vsi, wsi, Ri, wghti;
      double   s0m, s0p, s1m, s1p, s2m, s2p, t00, t10, t20, t0m, t0p, t1p;
      double   r00, r10, r20, r0p, r0m, r1p, f0, f1, f2, B0, B1, B2 ;
      double   IS, valr, vali ;
    
      switch ((int)Fct->Para[0]) {
    
      case _3D :
        
        switch (Element->ElementSource->Type) {
          
        case TRIANGLE :
        case QUADRANGLE :  
          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] ;
          
          valr = 0. ;
          vali = 0. ;
          IS = 0. ;
    
          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[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];
    	
    	l2[0] = (xs[2]-xs[0])*u[0] + (ys[2]-ys[0])*u[1] + (zs[2]-zs[0])*u[2]; /*u2 coordinate*/
    	l2[1] = (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-l2[0]) + yl*l2[1] ) / b;   s0p = s0m + b;
    	s1p = ( xl * l2[0] + yl * l2[1] ) / c;   s1m = s1p - c; s2m = - xl; s2p = a - xl;
    	
    	/*distance observation point projection on triangle plane to triangle local vertices*/
    	t00 = (yl * (l2[0]-a) + l2[1] * (a-xl)) / b;
    	t10 = (xl * l2[1] - yl * l2[0]) / c;
    	t20 = yl;
    	t0m = sqrt((a-xl)*(a-xl) + yl*yl); 
    	t0p = sqrt((l2[0]-xl)*(l2[0]-xl) + (l2[1]-yl)*(l2[1]-yl)); 
    	t1p = sqrt(xl*xl + yl*yl);
    	
    	/* minimum distances from the observation point to each triangle side*/
    	r00 = sqrt(SQU(t00) + SQU(zl));
    	r10 = sqrt(SQU(t10) + SQU(zl));
    	r20 = sqrt(SQU(t20) + SQU(zl));
    	
    	/* distances from observation point to the vertices*/
    	r0p = sqrt(SQU(t0p) + SQU(zl));      
    	r0m = sqrt(SQU(t0m) + SQU(zl));
    	r1p = sqrt(SQU(t1p) + SQU(zl));
    	
    	/* intermediate functions */
            
    	f0 = r00 <= EPSILON ? 0 : log((r0p + s0p) / (r0m + s0m));
    	f1 = r10 <= EPSILON ? 0 : log((r1p + s1p) / (r0p + s1m));
    	f2 = r20 <= EPSILON ? 0 : log((r0m + s2p) / (r1p + s2m));
    	
    	B0 = fabs(t00) <= EPSILON ? 0 : atan(t00*s0p/(SQU(r00)+fabs(zl)*r0p))-atan(t00*s0m/(SQU(r00)+fabs(zl)*r0m));
    	B1 = fabs(t10) <= EPSILON ? 0 : atan(t10*s1p/(SQU(r10)+fabs(zl)*r1p))-atan(t10*s1m/(SQU(r10)+fabs(zl)*r0p));
    	B2 = fabs(t20) <= EPSILON ? 0 : atan(t20*s2p/(SQU(r20)+fabs(zl)*r0m))-atan(t20*s2m/(SQU(r20)+fabs(zl)*r1p));
    	d = a * l2[1]; /* Double aire a cause de normalization */
    	
    	IS +=  ONE_OVER_TWO_PI * (-fabs(zl)*(B0+B1+B2) + t00*f0+t10*f1+t20*f2)/d; /* 1/r integral solution*/
    	
    	/* Gauss Numerical Integration of (exp(Fct->Para[1]*R)-1)/R */
    	for (i_IntPoint = 1; i_IntPoint <= NGT_Points; i_IntPoint++){
    	  Gauss_Triangle(NGT_Points,i_IntPoint,&us,&vs,&ws,&wghti);
    	  usi = u[0]*us + u[1]*vs + u[2]*ws ;
    	  vsi = v[0]*us + v[1]*vs + v[2]*ws ;
    	  wsi = n[0]*us + n[1]*vs + n[2]*ws ;
    	  Ri = sqrt( SQU(xl-usi) + SQU(yl-vsi) + SQU(zl-wsi) ) ;
    	  
    	  valr += Ri > EPSILON ? (wghti*(cos(Fct->Para[1]*Ri)-1)/Ri): 0 ;
    	  vali += Ri > EPSILON ? (-wghti*sin(Fct->Para[1]*Ri)/Ri): (-wghti * Fct->Para[1]); 
    	}
    	
    	valr = d * valr/2;
    	vali = d * vali/2;
    	
    	Val->Val[0] = IS + valr ;
    	Val->Val[MAX_DIM] = vali ; /* Imaginary part. Numerical integral */
    	
    	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];
    	}
          }
          if (j == 0){ Val->Val[0] = (Val->Val[0])/2; }
          Val->Type = SCALAR;
          break ;
          
        default :
          Message::Error("Unknown Element Type (%s) for 'GF_HelmholtzxForm'",
    		 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
        }
        break ;
        
      default :
        Message::Error("Unknown Dimension (%d) for 'GF_HelmholtzxForm'", 
    	       (int)Fct->Para[0]);
        
      }
    }