Commit 62aaf4af by Christophe Geuzaine

Merge branch 'minorFix' into 'master'

Minor fix (pp analytics + LIN_2/TRI_2/TET_2) See merge request !42
parents c91d487c cbdb31e5
Pipeline #1848 passed with stage
in 9 minutes 49 seconds
......@@ -1143,7 +1143,7 @@ struct StringXFunction2Nbr F_Function[] = { /* #Par #Arg */
{"AcousticFieldSoftCylinderABC", (CAST)F_AcousticFieldSoftCylinderABC, 5, 1 },
{"DrAcousticFieldSoftCylinder", (CAST)F_DrAcousticFieldSoftCylinder, 2, 1 },
{"RCSSoftCylinder", (CAST)F_RCSSoftCylinder, 2, 1 },
{"AcousticFieldHardCylinder", (CAST)F_AcousticFieldHardCylinder, 2, 1 },
{"AcousticFieldHardCylinder", (CAST)F_AcousticFieldHardCylinder, -1, 1 },
{"AcousticFieldHardCylinderABC", (CAST)F_AcousticFieldHardCylinderABC, 5, 1 },
{"DthetaAcousticFieldHardCylinder", (CAST)F_DthetaAcousticFieldHardCylinder, 2, 1 },
{"RCSHardCylinder", (CAST)F_RCSHardCylinder, 2, 1 },
......
......@@ -1686,7 +1686,7 @@ void F_AcousticFieldSoftCylinder(F_ARG)
std::complex<double> I(0,1);
double vr=0, vi=0;
#if defined(_OPENMP)
//#pragma omp parallel for reduction(+: vr,vi)
#pragma omp parallel for reduction(+: vr,vi)
#endif
for(int n = nStart ; n < nEnd ; n++){
std::complex<double> HnkR( jn(n,kR), yn(n,kR) );
......@@ -1944,49 +1944,48 @@ void F_RCSSoftCylinder(F_ARG)
void F_AcousticFieldHardCylinder(F_ARG)
{
cplx I = {0.,1.}, Hnkr, dHnkR, tmp, *HnkRtab;
double k, R, r, kr, kR, theta, cost ;
int n, ns ;
theta = atan2(A->Val[1], A->Val[0]) ;
r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]) ;
k = Fct->Para[0] ;
R = Fct->Para[1] ;
kr = k*r;
kR = k*R;
V->Val[0] = 0.;
V->Val[MAX_DIM] = 0. ;
ns = (int)k + 10;
HnkRtab = (cplx*)Malloc(ns*sizeof(cplx));
double theta = atan2(A->Val[1], A->Val[0]);
double r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]);
double k = Fct->Para[0];
double R = Fct->Para[1];
double kr = k*r;
double kR = k*R;
for (n = 0 ; n < ns ; n++){
HnkRtab[n].r = jn(n,kR);
HnkRtab[n].i = yn(n,kR);
// 3rd parameter : change the incidence angle
if(Fct->NbrParameters > 2){
double thetaInc = Fct->Para[2];
theta += thetaInc;
}
for (n = 0 ; n < ns ; n++){
Hnkr.r = jn(n,kr);
Hnkr.i = yn(n,kr);
dHnkR = DHn(HnkRtab, n, kR);
tmp = Cdiv( Cprod( Cpow(I,n) , Cprodr( dHnkR.r, Hnkr) ) , dHnkR );
cost = cos(n*theta);
V->Val[0] += cost * tmp.r * (!n ? 0.5 : 1.);
V->Val[MAX_DIM] += cost * tmp.i * (!n ? 0.5 : 1.);
// 4th/5th parameters : change the range of modes
int nStart = 0;
int nEnd = (int)(kR) + 10;
if(Fct->NbrParameters > 3){
int nStartNew = Fct->Para[3];
int nEndNew = Fct->Para[4];
nStart = nStartNew;
nEnd = nEndNew;
}
std::complex<double> *HnkRtab;
HnkRtab = new std::complex<double>[nEnd];
for (int n=nStart; n<nEnd; n++){
HnkRtab[n] = std::complex<double>(jn(n,kR),yn(n,kR));
}
Free(HnkRtab);
std::complex<double> I(0,1), val;
for (int n=nStart; n<nEnd; n++){
std::complex<double> Hnkr( jn(n,kr), yn(n,kr) );
std::complex<double> dHnkR = (!n ? -HnkRtab[1] : HnkRtab[n-1] - (double)n/kR * HnkRtab[n]);
std::complex<double> tmp1 = std::pow(I,n) * Hnkr/dHnkR;
double tmp2 = - (!n ? 1. : 2.) * cos(n*theta) * std::real(dHnkR);
val += tmp1 * tmp2;
}
V->Val[0] *= -2;
V->Val[MAX_DIM] *= -2;
delete HnkRtab;
V->Val[0] = std::real(val);
V->Val[MAX_DIM] = std::imag(val);
V->Type = SCALAR ;
}
......
......@@ -103,6 +103,7 @@ void Geo_CreateNormal(int Type, double *x, double *y, double *z, double *N)
switch (Type) {
case LINE :
case LINE_2 :
nx = y[1] - y[0] ;
ny = x[0] - x[1] ;
norm = sqrt(SQU(nx)+SQU(ny)) ;
......@@ -112,7 +113,9 @@ void Geo_CreateNormal(int Type, double *x, double *y, double *z, double *N)
break ;
case TRIANGLE :
case TRIANGLE_2 :
case QUADRANGLE :
case QUADRANGLE_2 :
x1x0 = x[1] - x[0] ;
y1y0 = y[1] - y[0] ;
z1z0 = z[1] - z[0] ;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment