#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 <getdp@geuz.org>.
*
* 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 ;
tx1 = xxs * c1 ;
ty1 = yys * c1 ;
tx2 = r_2 * xxs/rxy_2 ;
ty2 = r_2 * yys/rxy_2 ;
r_l_2 = -1/r/r;
for (l = 0 ; l < Nd ; l++){
for (m = -l ; m <= l ; m++){
Plm = Legendre(l, m, cTheta) ;
dPlm = dLegendre(l, m, cTheta) ;
sPhi = sin(m*Phi) ;
cPhi = cos(m*Phi) ;
cr = r_l_2/Factorial(l+m) ;
V[l*(l+1)+m].Type = VECTOR ;
V[l*(l+1)+m].Val[0] = cr *( Plm * m * ty2 * sPhi-
dPlm * tx1 * cPhi + l * xxs * Plm * cPhi) ;
V[l*(l+1)+m].Val[MAX_DIM ] = -cr *(-Plm * m * ty2 * cPhi-
dPlm * tx1 * sPhi + l * xxs * Plm * sPhi) ;
V[l*(l+1)+m].Val[1] = cr *(-Plm * m * tx2 * sPhi -
dPlm * ty1 * cPhi + l * yys * Plm * cPhi) ;
V[l*(l+1)+m].Val[MAX_DIM+1] = -cr *( Plm * m * tx2 * cPhi-
dPlm * ty1 * sPhi + l * yys * Plm * sPhi) ;
V[l*(l+1)+m].Val[2] = cr *(dPlm * srxy * cPhi + l * zzs * Plm * cPhi) ;
V[l*(l+1)+m].Val[MAX_DIM+2] = -cr *(dPlm * srxy * sPhi + l * zzs * Plm * sPhi) ;
}
r_l_2 *= 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)) ;
if(!r) Msg(GERROR, "1/r in 'GF_GradLaplace (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(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_GradLaplace' (%d)", Current.FMM.Flag_GF);
break;
}
break;
default :
Msg(GERROR, "Bad Parameter for 'GF_GradLaplace' (%d)", (int)Fct->Para[0]);
break;
}
V->Type = VECTOR ;
V->Val[MAX_DIM+0] = 0. ;
V->Val[MAX_DIM+1] = 0. ;
V->Val[MAX_DIM+2] = 0. ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G F _ N P x G r a d L a p l a c e */
/* ------------------------------------------------------------------------ */
void GF_NPxGradLaplace (F_ARG) {
int i, Nd, k, l, m ;
double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0, xxs, yys, zzs, a, b, c, n ;
double r, cte, cr, phi_r, sphi_r, cphi_r, r_l_2, r_l, r_l_1 ;
double Theta, Phi, sTheta, cTheta, Plm, dPlm, sPhi, cPhi ;
double r_2, rxy, rxy_2, srxy, c1, tx1, ty1, tx2, ty2 ;
GetDP_Begin("GF_NPxGradLaplace") ; /* It computes the Scalar product Normal[] * GradLaplace */
V->Type = SCALAR ;
V->Val[MAX_DIM] = 0. ;
if ((Current.Element->Num == Current.ElementSource->Num) && !Flag_FMM) {
V->Val[0 ] = 0. ;
V->Val[MAX_DIM] = 0. ;
GetDP_End ;
}
switch((int)Fct->Para[0]){
case _2D :
switch(Current.FMM.Flag_GF){
case FMM_DIRECT : /* Normal on Element */
x1x0 = Current.Element->x[1] - Current.Element->x[0] ;
y1y0 = Current.Element->y[1] - Current.Element->y[0] ;
xxs = Current.x-Current.xs ;
yys = Current.y-Current.ys ;
r = SQU(xxs)+SQU(yys) ;
if(!r) V->Val[0] = 0 ;
else
V->Val[0] = - ONE_OVER_TWO_PI * (y1y0 * xxs - x1x0 * yys)
/ (sqrt(SQU(x1x0)+SQU(y1y0)) * r) ;
V->Val[MAX_DIM] = 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 ;
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 :/* 2D */
Nd = Current.FMM.NbrDir ;
x1x0 = Current.Element->x[1] - Current.Element->x[0] ;
y1y0 = Current.Element->y[1] - Current.Element->y[0] ;
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_NPxGradLaplace (Disaggregation - 2D)' Integration point == Group center") ;
cte = -1/r/r/sqrt(SQU(x1x0)+SQU(y1y0)) ;
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 = SCALAR ;
V[i].Val[0] = ( y1y0 * (xxs * cphi_r + yys * sphi_r) - x1x0 * (yys * cphi_r - xxs * sphi_r)) * cte ;
V[i].Val[MAX_DIM] = (-y1y0 * (yys * cphi_r - xxs * sphi_r) - x1x0 * (xxs * cphi_r + yys * sphi_r)) * cte ;
}
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_NPxGradLaplace' (%d)", Current.FMM.Flag_GF);
break;
}
break ;
case _3D :
switch(Current.FMM.Flag_GF){
case FMM_DIRECT : /* Normal on Element */
x1x0 = Current.Element->x[1] - Current.Element->x[0] ;
y1y0 = Current.Element->y[1] - Current.Element->y[0] ;
z1z0 = Current.Element->z[1] - Current.Element->z[0] ;
x2x0 = Current.Element->x[2] - Current.Element->x[0] ;
y2y0 = Current.Element->y[2] - Current.Element->y[0] ;
z2z0 = Current.Element->z[2] - Current.Element->z[0] ;
a = y1y0 * z2z0 - z1z0 * y2y0 ;
b = z1z0 * x2x0 - x1x0 * z2z0 ;
c = x1x0 * y2y0 - y1y0 * x2x0 ;
xxs = Current.x-Current.xs ;
yys = Current.y-Current.ys ;
zzs = Current.z-Current.zs ;
V->Val[0] = - ONE_OVER_FOUR_PI * (a * xxs + b * yys + c * zzs)
/ ( (sqrt(SQU(a)+SQU(b)+SQU(c)) *
CUB(sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)))) ) ;
V->Val[MAX_DIM] = 0. ;
break ;
case FMM_AGGREGATION : /* SCALAR Grad 3D -> VECTOR in Disaggregation */
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 : /* Grad 3D -> VECTOR */
x1x0 = Current.Element->x[1] - Current.Element->x[0] ;
y1y0 = Current.Element->y[1] - Current.Element->y[0] ;
z1z0 = Current.Element->z[1] - Current.Element->z[0] ;
x2x0 = Current.Element->x[2] - Current.Element->x[0] ;
y2y0 = Current.Element->y[2] - Current.Element->y[0] ;
z2z0 = Current.Element->z[2] - Current.Element->z[0] ;
a = y1y0 * z2z0 - z1z0 * y2y0 ; /* Normal */
b = z1z0 * x2x0 - x1x0 * z2z0 ;
c = x1x0 * y2y0 - y1y0 * x2x0 ;
n = sqrt(SQU(a)+SQU(b)+SQU(c)) ;
a /= n ;/* Normalized Normal on Element */
b /= n ;
c /= n ;
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_NPxGradLaplace (Disaggregation - 3D)' r %e ",r) ;
rxy_2 = SQU(xxs)+SQU(yys) ;
rxy = sqrt(rxy_2) ;
if(!rxy) Msg(GERROR, "1/0 in 'GF_NPxGradLaplace (Disaggregation - 3D)' rxy %e",rxy ) ;
Theta = atan2(rxy, zzs) ;
Phi = atan2(yys, xxs) ;
cTheta = cos(Theta) ;
sTheta = sin(Theta) ;
srxy = sTheta * rxy ;
c1 = sTheta * zzs/rxy ;
tx1 = xxs * c1 ;
ty1 = yys * c1 ;
tx2 = r_2 * xxs/rxy_2 ;
ty2 = r_2 * yys/rxy_2 ;
r_l_2 = -1/r/r;
for (l = 0 ; l < Nd ; l++){
for (m = -l ; m <= l ; m++){
Plm = Legendre(l, m, cTheta) ;
dPlm = dLegendre(l, m, cTheta) ;
sPhi = sin(m*Phi) ;
cPhi = cos(m*Phi) ;
cr = r_l_2/Factorial(l+m) ;
V[l*(l+1)+m].Type = SCALAR ;
V[l*(l+1)+m].Val[0] = a * cr *( Plm * m * ty2 * sPhi - dPlm * tx1 * cPhi + l * xxs * Plm * cPhi) +
b * cr *(-Plm * m * tx2 * sPhi - dPlm * ty1 * cPhi + l * yys * Plm * cPhi) +
c * cr *(dPlm * srxy * cPhi + l * zzs * Plm * cPhi) ;
V[l*(l+1)+m].Val[MAX_DIM] = -a * cr *(-Plm * m * ty2 * cPhi - dPlm * tx1 * sPhi + l * xxs * Plm * sPhi)-
b * cr *( Plm * m * tx2 * cPhi - dPlm * ty1 * sPhi + l * yys * Plm * sPhi)-
c * cr *(dPlm * srxy * sPhi + l * zzs * Plm * sPhi) ;
}
r_l_2 *= r ;
}
break ;
case FMM_TRANSLATION :
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)) ;
if(!r) Msg(GERROR, "1/0 in 'GF_NPxGradLaplace (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_NPxGradLaplace' (%d)", Current.FMM.Flag_GF);
break;
}
break;
default :
Msg(GERROR, "Bad Parameter for 'GF_NPxGradLaplace' (%d)", (int)Fct->Para[0]);
break;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G F _ N S x G r a d L a p l a c e */
/* ------------------------------------------------------------------------ */
void GF_NSxGradLaplace (F_ARG) {
double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0, xxs, yys, zzs, a, b, c ;
GetDP_Begin("GF_NSxGradLaplace");
V->Type = SCALAR ;
V->Val[MAX_DIM] = 0. ;
if (Current.Element->Num == Current.ElementSource->Num) {
V->Val[0] = 0. ;
V->Val[MAX_DIM] = 0. ;
GetDP_End ;
}
switch((int)Fct->Para[0]){
case _2D :
x1x0 = Current.ElementSource->x[1] - Current.ElementSource->x[0] ;
y1y0 = Current.ElementSource->y[1] - Current.ElementSource->y[0] ;
xxs = Current.x-Current.xs ;
yys = Current.y-Current.ys ;
V->Val[0] = ONE_OVER_TWO_PI * (y1y0 * xxs - x1x0 * yys)
/ (sqrt(SQU(x1x0)+SQU(y1y0)) * (SQU(xxs)+SQU(yys)) ) ;
V->Val[MAX_DIM] = 0. ;
break ;
case _3D :
x1x0 = Current.ElementSource->x[1] - Current.ElementSource->x[0] ;
y1y0 = Current.ElementSource->y[1] - Current.ElementSource->y[0] ;
z1z0 = Current.ElementSource->z[1] - Current.ElementSource->z[0] ;
x2x0 = Current.ElementSource->x[2] - Current.ElementSource->x[0] ;
y2y0 = Current.ElementSource->y[2] - Current.ElementSource->y[0] ;
z2z0 = Current.ElementSource->z[2] - Current.ElementSource->z[0] ;
a = y1y0 * z2z0 - z1z0 * y2y0 ;
b = z1z0 * x2x0 - x1x0 * z2z0 ;
c = x1x0 * y2y0 - y1y0 * x2x0 ;
xxs = Current.x-Current.xs ;
yys = Current.y-Current.ys ;
zzs = Current.z-Current.zs ;
V->Val[0] = ONE_OVER_FOUR_PI * (a * xxs + b * yys + c * zzs)
/ ( (sqrt(SQU(a)+SQU(b)+SQU(c)) *
CUB(sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)))) ) ;
V->Val[MAX_DIM] = 0. ;
break ;
default :
Msg(GERROR, "Bad Parameter for 'GF_NSxGradLaplace' (%d)", (int)Fct->Para[0]);
break;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G F _ A p p r o x i m a t e L a p l a c e */
/* ------------------------------------------------------------------------ */
void GF_ApproximateLaplace (F_ARG) {
GetDP_Begin("GF_ApproxilateLaplace");
Msg(GERROR, "The Approximate Integral Kernels can only be Integrated Analytically");
GetDP_End ;
}
#undef F_ARG
syntax highlighted by Code2HTML, v. 0.9.1