#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 <getdp@geuz.org>.
 *
 * 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 :
    
    switch (Element->ElementSource->Type) {
      
    case POINT :
      Val->Type = VECTOR ;
      
      if (Element->Num == Element->ElementSource->Num) {
	Val->Val[0] = Val->Val[1] = Val->Val[2] = 0. ;
	GetDP_End ;
      }
      
      xxs = x - Element->ElementSource->x[0] ;
      yys = y - Element->ElementSource->y[0] ;
      r2 = SQU(xxs)+SQU(yys) ;
      
      if (r2 > EPSILON2) {
	Val->Val[0] = - ONE_OVER_TWO_PI * xxs / r2 ;
	Val->Val[1] = - ONE_OVER_TWO_PI * yys / r2 ;
	Val->Val[2] = 0. ;
      }
      else {
	Val->Val[0] = Val->Val[1] = Val->Val[2] = 0. ;
      }
      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) {
	cglt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	Val->Type = VECTOR ;
	Val->Val[0] = valr ;
	Val->Val[1] = vali ;
	Val->Val[2] = 0. ;
      }
      if(xFunctionBF == (CAST)BF_GradNode) {
	DetJ = (xs[2]-xs[0]) * (ys[1]-ys[0]) - (xs[1]-xs[0]) * (ys[2]-ys[0]) ;
	switch(NumEntity){
	case 1 : mx = (ys[2]-ys[1])/DetJ ; my = (xs[1]-xs[2])/DetJ ; break;
	case 2 : mx = (ys[0]-ys[2])/DetJ ; my = (xs[2]-xs[0])/DetJ ; break;
	case 3 : mx = (ys[1]-ys[0])/DetJ ; my = (xs[0]-xs[1])/DetJ ; break;
	}
	cglt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	Val->Type = SCALAR ;
	Val->Val[0] = my*valr - mx*vali ; 
      }
      else{
	Msg(GERROR, "Unknown Basis Function Type for 'GF_GradLaplacexForm'");
      }
      break ;

    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_GradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));
    }
    break ;
    
    
  case _3D :
    
    switch (Element->ElementSource->Type) {
      
    case LINE :
      Val->Type = VECTOR ;

      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) ;
      
      I1 = 2./(4.*a*c-b*b) * ( (2.*a+b)/sqrt(a+b+c) - b/sqrt(c) ) ;
      I2 = 2./(-4.*a*c+b*b) * ( (2.*c+b)/sqrt(a+b+c) - 2.*sqrt(c) ) ;
      a2 = sqrt(a) ;
      
      Val->Val[0] = ONE_OVER_FOUR_PI * ( (xs[0]-x) * I1 + (xs[1]-xs[0]) * I2 ) * a2 ;
      Val->Val[1] = ONE_OVER_FOUR_PI * ( (ys[0]-y) * I1 + (ys[1]-ys[0]) * I2 ) * a2 ;
      Val->Val[2] = ONE_OVER_FOUR_PI * ( (zs[0]-z) * I1 + (zs[1]-zs[0]) * I2 ) * a2 ;
      break ;

    case TRIANGLE : 
      Val->Type = VECTOR ;

      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(xFunctionBF == (CAST)BF_Volume) Type_Int = 1 ;
      if(xFunctionBF == (CAST)BF_Node)   Type_Int = 2 ;

      /* 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];
	
      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]);

      area = a * v2/2 ;/* Triangle area */
      if (!zl) zl = sqrt(area) * 1e-15 ;

      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*/
      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 ;
     
      r00 = sqrt(r00_2);
      r10 = sqrt(r10_2);
      r20 = sqrt(r20_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);
      
      EPS = EPSILON*(fabs(s0m)+fabs(s0p));

      B0 = (r00 <= EPS) ? 0. : atan(t00*s0p/(r00_2+fabs(zl)*r0p))-atan(t00*s0m/(r00_2+fabs(zl)*r0m));
      f20 = ((r0m + s0m) <= EPS) ? log(s0m/s0p) : log((r0p + s0p) / (r0m + s0m)) ;

      EPS = EPSILON*(fabs(s1m)+fabs(s1p)) ; 
      B1 = (r10 <=EPS) ? 0. : atan(t10*s1p/(r10_2+fabs(zl)*r1p))-atan(t10*s1m/(r10_2+fabs(zl)*r0p));
      f21 = ((r0p + s1m)<=EPS) ? log(s1m/s1p) : log((r1p + s1p) / (r0p + s1m));

      EPS = EPSILON*(fabs(s2m)+fabs(s2p)) ; 
      B2 =  (r20 <= EPS) ? 0. : atan(t20*s2p/(r20_2+fabs(zl)*r0m))-atan(t20*s2m/(r20_2+fabs(zl)*r1p));
      f22 = ((r1p + s2m)< EPS) ? log(s2m/s2p): log((r0m + s2p) / (r1p + s2m));
      
      B = B0 + B1 + B2 ; 

      s0[0] = (xs[2] - xs[1])/b ;
      s0[1] = (ys[2] - ys[1])/b ;
      s0[2] = (zs[2] - zs[1])/b ;

      s1[0] = (xs[0] - xs[2])/c ;
      s1[1] = (ys[0] - ys[2])/c ;
      s1[2] = (zs[0] - zs[2])/c ;

      m0[0] = s0[1] * n[2] - s0[2]* n[1] ;
      m0[1] = s0[2] * n[0] - s0[0]* n[2] ;
      m0[2] = s0[0] * n[1] - s0[1]* n[0] ;

      m1[0] = s1[1]* n[2] - s1[2] * n[1] ;
      m1[1] = s1[2]* n[0] - s1[0] * n[2] ;
      m1[2] = s1[0]* n[1] - s1[1] * n[0] ;
      
      m2[0] = u[1] * n[2] - u[2]* n[1] ;
      m2[1] = u[2] * n[0] - u[0]* n[2] ;
      m2[2] = u[0] * n[1] - u[1]* n[0] ;

      /* Grad(1/r) integral solution*/
      I[0] = -n[0] * THESIGN(zl) * B - (m0[0]*f20 + m1[0]*f21 + m2[0]*f22) ; 
      I[1] = -n[1] * THESIGN(zl) * B - (m0[1]*f20 + m1[1]*f21 + m2[1]*f22) ;
      I[2] = -n[2] * THESIGN(zl) * B - (m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;

      switch ( Type_Int ){
      case 1 : /* BF_Volume */	
	Val->Val[0] = I[0]/area ; 	
	Val->Val[1] = I[1]/area ; 	
	Val->Val[2] = I[2]/area ; 	
	break;
    
      case 2 : /* BF_Node */
	if (!v2 ) 	
	  Msg(GERROR,"1/0 in GF_LaplacexForm (case _3D TRIANGLE) v2 %e", v2);

	f0[0] = s0[0] * t00 * f20 - m0[0]*(r0p-r0m) ; /* fi */
	f0[1] = s0[1] * t00 * f20 - m0[1]*(r0p-r0m) ;
	f0[2] = s0[2] * t00 * f20 - m0[2]*(r0p-r0m) ;

	f1[0] = s1[0] * t10 * f21 - m1[0]*(r1p-r0p) ;
	f1[1] = s1[1] * t10 * f21 - m1[1]*(r1p-r0p) ;
	f1[2] = s1[2] * t10 * f21 - m1[2]*(r1p-r0p) ;

	f2[0] = u[0] * t20 * f22 - m2[0]*(r0m-r1p) ;
	f2[1] = u[1] * t20 * f22 - m2[1]*(r0m-r1p) ;
	f2[2] = u[2] * t20 * f22 - m2[2]*(r0m-r1p) ;  

	umf2i = u[0]*(m0[0]*f20 + m1[0]*f21 + m2[0]*f22) +
	        u[1]*(m0[1]*f20 + m1[1]*f21 + m2[1]*f22) + 
	        u[2]*(m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;

	us0 = u[0] * s0[0] + u[1] * s0[1] + u[2] * s0[2] ; 
	us1 = u[0] * s1[0] + u[1] * s1[1] + u[2] * s1[2] ; 
	us2 = u[0] *  u[0] + u[1] *  u[1] + u[2] *  u[2] ; 

	vmf2i = v[0]*(m0[0]*f20 + m1[0]*f21 + m2[0]*f22) + 
	        v[1]*(m0[1]*f20 + m1[1]*f21 + m2[1]*f22) + 
	        v[2]*(m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;

	vs0 = v[0] * s0[0] + v[1] * s0[1] + v[2] * s0[2] ; 
	vs1 = v[0] * s1[0] + v[1] * s1[1] + v[2] * s1[2] ; 
	vs2 = v[0] *  u[0] + v[1] *  u[1] + v[2] *  u[2] ; 

	B *= fabs(zl);
	umf2i *= zl ;
	vmf2i *= zl ;

	Iua[0] = n[0] * umf2i - B * u[0] + f0[0] * us0 + f1[0] * us1 + f2[0] * us2 ;
	Iua[1] = n[1] * umf2i - B * u[1] + f0[1] * us0 + f1[1] * us1 + f2[1] * us2 ;
	Iua[2] = n[2] * umf2i - B * u[2] + f0[2] * us0 + f1[2] * us1 + f2[2] * us2 ;
	
	Iva[0] = n[0] * vmf2i - B * v[0] + f0[0] * vs0 + f1[0] * vs1 + f2[0] * vs2 ;
	Iva[1] = n[1] * vmf2i - B * v[1] + f0[1] * vs0 + f1[1] * vs1 + f2[1] * vs2 ;
	Iva[2] = n[2] * vmf2i - B * v[2] + f0[2] * vs0 + f1[2] * vs1 + f2[2] * vs2 ;
	
	switch(NumEntity){
	case 1 :  
	  N10 = 1 - xl/a + (u2/a -1) * yl/v2 ;
	  Val->Val[0] = N10 * I[0]  - Iua[0]/a + (u2/a-1) * Iva[0]/v2 ;
	  Val->Val[1] = N10 * I[1]  - Iua[1]/a + (u2/a-1) * Iva[1]/v2 ;
	  Val->Val[2] = N10 * I[2]  - Iua[2]/a + (u2/a-1) * Iva[2]/v2 ;
	  break;
	case 2 : 
	  N20 = xl/a - u2/a * yl/v2 ;
	  Val->Val[0] = N20 * I[0] + Iua[0]/a - u2/a * Iva[0]/v2 ;
	  Val->Val[1] = N20 * I[1] + Iua[1]/a - u2/a * Iva[1]/v2 ;
	  Val->Val[2] = N20 * I[2] + Iua[2]/a - u2/a * Iva[2]/v2 ;
	  break;
	case 3 :  
	  N30 = yl/v2 ;
	  Val->Val[0] = N30 * I[0] + Iva[0]/v2 ;
	  Val->Val[1] = N30 * I[1] + Iva[1]/v2 ;
	  Val->Val[2] = N30 * I[2] + Iva[2]/v2 ;
	  break;
	}
	break;
      }
      
      Val->Val[0] *= ONE_OVER_FOUR_PI ;
      Val->Val[1] *= ONE_OVER_FOUR_PI ;
      Val->Val[2] *= ONE_OVER_FOUR_PI ;      
      break ;
      
    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_GradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));
    }
    break ;
    
  default :
    Msg(GERROR, "Unknown Dimension (%d) for 'GF_GradLaplacexForm'",
	(int)Fct->Para[0]);
    
  }

  GetDP_End ;
}


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

void GF_NPxGradLaplacexForm (F_ARG2) {

  double  xs[MAX_NODES], ys[MAX_NODES] ;
  double  xp[MAX_NODES], yp[MAX_NODES], N[3] ;
  int     Type_Int;
  double  a, b, c, d, m, n, Jp, i1, Is, I1=0 ;
  struct Value ValGrad ;

  GetDP_Begin("GF_NPxGradLaplacexForm");

  Val->Type = SCALAR ;
  Val->Val[MAX_DIM] = 0.0 ;
  
  if (Element->Num == Element->ElementSource->Num) {
    Val->Val[0] = 0.0 ;
    GetDP_End ;
  }
  
  switch ((int)Fct->Para[0]) {
    
  case _2D : 

    switch (Element->ElementSource->Type) {

    case LINE :	  
      if (Element->Type != LINE)
	Msg(GERROR, "GF_NPxGradLaplacexForm not ready for mixed geometrical elements");

      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) {
	if ((x == xs[0]) && (y == ys[0]))
	  Type_Int = 1 ; 
	else if ((x == xs[1]) && (y == ys[1]))
	  Type_Int = 2 ; 
	else                                    
	  Type_Int = 3 ; 
	
	xp[0] = Element->x[0] ; yp[0] = Element->y[0] ;
	xp[1] = Element->x[1] ; yp[1] = Element->y[1] ;
	
	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 = 4.*a*c - b*b ;
	
	switch (Type_Int) {	  
	case 1 :  
	case 2 :  
	  Msg(GERROR, "Degenerate case not done in 'GF_NPxGradLaplacexForm'");
	  break ;	
	case 3 :  
	  if (fabs(d) < EPSILON2) {
	    I1 = 0.0 ;
	  }
	  else {
	    if(d<0) Msg(GERROR, "Unexpected value in 'GF_NPxGradLaplacexForm'");
	    i1 = sqrt(d) ;
	    Is = 2. / i1 * (atan((2.*a+b)/i1) -  atan(b/i1)) ;
	    Jp = sqrt(SQU(xp[0]-xp[1])+SQU(yp[0]-yp[1])) ;
	    m = ((ys[0]-ys[1]) * (xp[0]-xp[1]) + (xs[0]-xs[1]) * (yp[1]-yp[0])) / Jp ;
	    n = ((yp[1]-yp[0]) * (x-xs[0]) + (xp[0]-xp[1]) * (y-ys[0])) / Jp ; 	
	    I1 = m /(2.*a) * log((a+b)/c+1.) + (n - m*b/(2.*a)) * Is ;	   
	  }
	  break ;
	}
	Val->Val[0] = - ONE_OVER_TWO_PI * I1 ; 
      }
      else {	
	Msg(GERROR, "Unknown Basis Function Type for 'GF_NPxGradLaplacexForm'");
      }
      break ;

    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_NPxGradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));     
    }
    break ;

  case _3D:
    switch (Element->ElementSource->Type) {
    case TRIANGLE :
      Geo_CreateNormal(Element->Type,
		       Element->x,Element->y,Element->z, N);

      GF_GradLaplacexForm(Element, Fct, xFunctionBF, NumEntity, x, y, z, &ValGrad) ;

      Val->Val[0] = N[0]*ValGrad.Val[0] + N[1]*ValGrad.Val[1] + N[2]*ValGrad.Val[2] ;
      break ;
    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_NPxGradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));  
    }
    break ;

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

  GetDP_End ;
}



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

void GF_NSxGradLaplacexForm (F_ARG2) {
    
  GetDP_Begin("GF_NSxGradLaplacexForm");

  Msg(GERROR, "Not done: 'GF_NSxGradLaplacexForm'");
  
  GetDP_End ;
}


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

void GF_ApproximateLaplacexForm (F_ARG2) {

  GetDP_Begin("GF_ApproxilateLaplacexForm");

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

  case 0 :
    GF_LaplacexForm(Element, Fct, (CAST)BF_Volume, 1, x, y, z, Val);
    break ;

  default :
    Msg(GERROR, "Bad Parameter Value in 'GF_ApproximateLaplacexForm'");
    break;

  }

  GetDP_End ;
}


#undef F_ARG2
#undef CAST