"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Functions/GF_Helmholtz.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

[ To the main getdp source changes report ]

GF_Helmholtz.cpp  (getdp-3.4.0-source.tgz):GF_Helmholtz.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
// //
// Contributor(s): // Contributor(s):
// Ruth Sabariego // Ruth Sabariego
// //
// g++ -std=c++11 on mingw does not define bessel functions // g++ -std=c++11 on mingw does not define bessel functions
#if defined (WIN32) && !defined(__CYGWIN__) #if defined(WIN32) && !defined(__CYGWIN__)
#undef __STRICT_ANSI__ #undef __STRICT_ANSI__
#endif #endif
#include <math.h> #include <math.h>
#include "GetDPConfig.h" #include "GetDPConfig.h"
#include "ProData.h" #include "ProData.h"
#include "GF.h" #include "GF.h"
#include "Cal_Value.h" #include "Cal_Value.h"
#include "Message.h" #include "Message.h"
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
#include "GeoData.h" #include "GeoData.h"
#endif #endif
#define SQU(a) ((a)*(a)) #define SQU(a) ((a) * (a))
#define CUB(a) ((a)*(a)*(a)) #define CUB(a) ((a) * (a) * (a))
#define ONE_OVER_FOUR_PI 7.9577471545947668E-02 #define ONE_OVER_FOUR_PI 7.9577471545947668E-02
extern struct CurrentData Current ; extern struct CurrentData Current;
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G F _ H e l m h o l t z */ /* G F _ H e l m h o l t z */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void GF_Helmholtz(GF_ARG) void GF_Helmholtz(GF_ARG)
{ {
double r, kr ; double r, kr;
if(Current.NbrHar != 2) if(Current.NbrHar != 2)
Message::Error("Wrong Number of Harmonics in 'GF_Helmholtz'"); Message::Error("Wrong Number of Harmonics in 'GF_Helmholtz'");
V->Type = SCALAR ; V->Type = SCALAR;
switch((int)Fct->Para[0]){ switch((int)Fct->Para[0]) {
case DIM_2D : case DIM_2D:
r = sqrt(SQU(Current.x-Current.xs)+ r = sqrt(SQU(Current.x - Current.xs) + SQU(Current.y - Current.ys));
SQU(Current.y-Current.ys) ) ; if(!r) Message::Error("1/0 in 'GF_Helmholtz'");
if(!r) Message::Error("1/0 in 'GF_Helmholtz'") ; kr = Fct->Para[1] * r;
kr = Fct->Para[1]*r;
V->Val[0] = -y0(kr)/4 ;
V->Val[MAX_DIM] = -j0(kr)/4 ;
break ;
case DIM_3D :
r = sqrt(SQU(Current.x-Current.xs)+
SQU(Current.y-Current.ys)+
SQU(Current.z-Current.zs)) ;
if(!r) Message::Error("1/0 in 'GF_Helmholtz'") ;
kr = Fct->Para[1]*r;
V->Val[0] = ONE_OVER_FOUR_PI * cos(kr) / r ;
V->Val[MAX_DIM] = -ONE_OVER_FOUR_PI * sin(kr) / r ;
break ;
default : V->Val[0] = -y0(kr) / 4;
V->Val[MAX_DIM] = -j0(kr) / 4;
break;
case DIM_3D:
r = sqrt(SQU(Current.x - Current.xs) + SQU(Current.y - Current.ys) +
SQU(Current.z - Current.zs));
if(!r) Message::Error("1/0 in 'GF_Helmholtz'");
kr = Fct->Para[1] * r;
V->Val[0] = ONE_OVER_FOUR_PI * cos(kr) / r;
V->Val[MAX_DIM] = -ONE_OVER_FOUR_PI * sin(kr) / r;
break;
default:
Message::Error("Bad Parameter for 'GF_Helmholtz' (%d)", (int)Fct->Para[0]); Message::Error("Bad Parameter for 'GF_Helmholtz' (%d)", (int)Fct->Para[0]);
break; break;
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G F _ H e l m h o l t z T h i n W i r e */ /* G F _ H e l m h o l t z T h i n W i r e */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void GF_HelmholtzThinWire(GF_ARG) void GF_HelmholtzThinWire(GF_ARG)
{ {
double a , r, kr ; double a, r, kr;
if(Current.NbrHar != 2) if(Current.NbrHar != 2)
Message::Error("Wrong Number of Harmonics in 'GF_HelmholtzThinWire'"); Message::Error("Wrong Number of Harmonics in 'GF_HelmholtzThinWire'");
V->Type = SCALAR ; V->Type = SCALAR;
switch((int)Fct->Para[0]) {
case DIM_2D:
a = Fct->Para[2];
r =
sqrt(SQU(Current.x - Current.xs) + SQU(Current.y - Current.ys) + SQU(a));
if(!r) Message::Error("1/0 in 'GF_HelmholtzThinWire'");
kr = Fct->Para[1] * r;
switch((int)Fct->Para[0]){ V->Val[0] = -y0(kr) / 4;
case DIM_2D : V->Val[MAX_DIM] = -j0(kr) / 4;
a = Fct->Para[2] ; break;
r = sqrt(SQU(Current.x-Current.xs)+
SQU(Current.y-Current.ys)+SQU(a)) ; case DIM_3D:
if(!r) Message::Error("1/0 in 'GF_HelmholtzThinWire'") ; a = Fct->Para[2];
kr = Fct->Para[1]*r;
r = sqrt(SQU(Current.x - Current.xs) + SQU(Current.y - Current.ys) +
V->Val[0] = -y0(kr)/4 ; SQU(Current.z - Current.zs) + SQU(a));
V->Val[MAX_DIM] = -j0(kr)/4 ; if(!r) Message::Error("1/0 in 'GF_HelmholtzThinWire'");
break ;
kr = Fct->Para[1] * r;
case DIM_3D : V->Val[0] = ONE_OVER_FOUR_PI * cos(kr) / r;
a = Fct->Para[2] ; V->Val[MAX_DIM] = -ONE_OVER_FOUR_PI * sin(kr) / r;
break;
r = sqrt(SQU(Current.x-Current.xs)+
SQU(Current.y-Current.ys)+
SQU(Current.z-Current.zs)+SQU(a)) ;
if(!r) Message::Error("1/0 in 'GF_HelmholtzThinWire'") ;
kr = Fct->Para[1]*r;
V->Val[0] = ONE_OVER_FOUR_PI * cos(kr) / r ;
V->Val[MAX_DIM] = -ONE_OVER_FOUR_PI * sin(kr) / r ;
break ;
default : default:
Message::Error("Bad Parameter for 'GF_HelmholtzThinWire' (%d)", (int)Fct->Pa Message::Error("Bad Parameter for 'GF_HelmholtzThinWire' (%d)",
ra[0]); (int)Fct->Para[0]);
break; break;
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G F _ G r a d H e l m h o l t z */ /* G F _ G r a d H e l m h o l t z */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* the gradient is taken relative to the destination point (x,y,z) */ /* the gradient is taken relative to the destination point (x,y,z) */
void GF_GradHelmholtz(GF_ARG) void GF_GradHelmholtz(GF_ARG)
{ {
double xxs, yys, zzs, r, kr, k0r ; double xxs, yys, zzs, r, kr, k0r;
double c1, c2, cr, ci ; double c1, c2, cr, ci;
if(Current.NbrHar != 2) if(Current.NbrHar != 2)
Message::Error("Wrong Number of Harmonics in 'GF_GradHelmholtz'"); Message::Error("Wrong Number of Harmonics in 'GF_GradHelmholtz'");
V->Type = VECTOR ; V->Type = VECTOR;
switch((int)Fct->Para[0]){ switch((int)Fct->Para[0]) {
case DIM_2D : case DIM_2D:
xxs = Current.x-Current.xs ; xxs = Current.x - Current.xs;
yys = Current.y-Current.ys ; yys = Current.y - Current.ys;
r = sqrt(SQU(xxs)+SQU(yys)) ; r = sqrt(SQU(xxs) + SQU(yys));
k0r = Fct->Para[1]*r; k0r = Fct->Para[1] * r;
if (!r) Cal_ZeroValue(V); if(!r)
Cal_ZeroValue(V);
else { else {
c1 = Fct->Para[1]/4/r ; c1 = Fct->Para[1] / 4 / r;
cr = c1 * y1(k0r); cr = c1 * y1(k0r);
ci = c1 * j1(k0r); ci = c1 * j1(k0r);
V->Val[0] = xxs * cr ; V->Val[MAX_DIM ] = xxs * ci ; V->Val[0] = xxs * cr;
V->Val[1] = yys * cr ; V->Val[MAX_DIM+1] = yys * ci ; V->Val[MAX_DIM] = xxs * ci;
V->Val[1] = yys * cr;
V->Val[MAX_DIM + 1] = yys * ci;
} }
break ; break;
case DIM_3D : case DIM_3D:
xxs = Current.x-Current.xs ; xxs = Current.x - Current.xs;
yys = Current.y-Current.ys ; yys = Current.y - Current.ys;
zzs = Current.z-Current.zs ; zzs = Current.z - Current.zs;
r = sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)) ; r = sqrt(SQU(xxs) + SQU(yys) + SQU(zzs));
kr = Fct->Para[1] * r ; kr = Fct->Para[1] * r;
if (!r) Cal_ZeroValue(V); if(!r)
Cal_ZeroValue(V);
else { else {
c1 = - ONE_OVER_FOUR_PI / CUB(r) ; c1 = -ONE_OVER_FOUR_PI / CUB(r);
c2 = ONE_OVER_FOUR_PI * Fct->Para[1] / SQU(r) ; c2 = ONE_OVER_FOUR_PI * Fct->Para[1] / SQU(r);
cr = c1 * cos(kr) - c2 * sin(kr) ; cr = c1 * cos(kr) - c2 * sin(kr);
ci = -c1 * sin(kr) - c2 * cos(kr) ; ci = -c1 * sin(kr) - c2 * cos(kr);
V->Val[0] = xxs * cr ; V->Val[MAX_DIM ] = xxs * ci ; V->Val[0] = xxs * cr;
V->Val[1] = yys * cr ; V->Val[MAX_DIM+1] = yys * ci ; V->Val[MAX_DIM] = xxs * ci;
V->Val[2] = zzs * cr ; V->Val[MAX_DIM+2] = zzs * ci ; V->Val[1] = yys * cr;
V->Val[MAX_DIM + 1] = yys * ci;
V->Val[2] = zzs * cr;
V->Val[MAX_DIM + 2] = zzs * ci;
} }
break ;
default :
Message::Error("Bad Parameter for 'GF_GradHelmholtz' (%d)", (int)Fct->Para[0
]);
break; break;
default:
Message::Error("Bad Parameter for 'GF_GradHelmholtz' (%d)",
(int)Fct->Para[0]);
break;
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G F _ N P x G r a d H e l m h o l t z */ /* G F _ N P x G r a d H e l m h o l t z */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void GF_NPxGradHelmholtz(GF_ARG) void GF_NPxGradHelmholtz(GF_ARG)
{ {
#if !defined(HAVE_KERNEL) #if !defined(HAVE_KERNEL)
Message::Error("GF_NPxGradHelmholtz requires Kernel"); Message::Error("GF_NPxGradHelmholtz requires Kernel");
#else #else
double N[3] ; double N[3];
struct Value ValGrad ; struct Value ValGrad;
/* Vectorial product N[] /\ Grad G */ /* Vectorial product N[] /\ Grad G */
if(Current.NbrHar != 2) if(Current.NbrHar != 2)
Message::Error("Wrong Number of Harmonics in 'GF_NPxGradHelmholtz'"); Message::Error("Wrong Number of Harmonics in 'GF_NPxGradHelmholtz'");
V->Type = VECTOR ; V->Type = VECTOR;
if (Current.Element->Num == Current.ElementSource->Num) { if(Current.Element->Num == Current.ElementSource->Num) {
Cal_ZeroValue(V); Cal_ZeroValue(V);
return ; return;
} }
switch((int)Fct->Para[0]){ switch((int)Fct->Para[0]) {
case DIM_3D : case DIM_3D:
Geo_CreateNormal(Current.Element->Type, Geo_CreateNormal(Current.Element->Type, Current.Element->x,
Current.Element->x,Current.Element->y,Current.Element->z, N) Current.Element->y, Current.Element->z, N);
;
GF_GradHelmholtz(Fct, &ValGrad, &ValGrad);
GF_GradHelmholtz(Fct, &ValGrad, &ValGrad) ;
V->Val[0] = N[1] * ValGrad.Val[2] - N[2] * ValGrad.Val[1];
V->Val[0] = N[1]*ValGrad.Val[2] - N[2]*ValGrad.Val[1]; V->Val[1] = -N[0] * ValGrad.Val[2] + N[2] * ValGrad.Val[0];
V->Val[1] =-N[0]*ValGrad.Val[2] + N[2]*ValGrad.Val[0]; V->Val[2] = N[0] * ValGrad.Val[1] - N[1] * ValGrad.Val[0];
V->Val[2] = N[0]*ValGrad.Val[1] - N[1]*ValGrad.Val[0]; V->Val[MAX_DIM] =
V->Val[MAX_DIM ] = N[1]*ValGrad.Val[MAX_DIM+2] - N[2]*ValGrad.Val[MAX_DIM+1 N[1] * ValGrad.Val[MAX_DIM + 2] - N[2] * ValGrad.Val[MAX_DIM + 1];
]; V->Val[MAX_DIM + 1] =
V->Val[MAX_DIM+1] =-N[0]*ValGrad.Val[MAX_DIM+2] + N[2]*ValGrad.Val[MAX_DIM]; -N[0] * ValGrad.Val[MAX_DIM + 2] + N[2] * ValGrad.Val[MAX_DIM];
V->Val[MAX_DIM+2] = N[0]*ValGrad.Val[MAX_DIM+1] - N[1]*ValGrad.Val[MAX_DIM]; V->Val[MAX_DIM + 2] =
break ; N[0] * ValGrad.Val[MAX_DIM + 1] - N[1] * ValGrad.Val[MAX_DIM];
break;
default : default:
Message::Error("Bad Parameter for 'GF_NPxGradHelmholtz' (%d)", (int)Fct->Par Message::Error("Bad Parameter for 'GF_NPxGradHelmholtz' (%d)",
a[0]); (int)Fct->Para[0]);
break; break;
} }
#endif #endif
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G F _ N S x G r a d H e l m h o l t z */ /* G F _ N S x G r a d H e l m h o l t z */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void GF_NSxGradHelmholtz(GF_ARG) void GF_NSxGradHelmholtz(GF_ARG)
{ {
double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0, xxs, yys, zzs, r ; double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0, xxs, yys, zzs, r;
double nx, ny, nz, n, c1, c2, cr, ci ; double nx, ny, nz, n, c1, c2, cr, ci;
if(Current.NbrHar != 2) if(Current.NbrHar != 2)
Message::Error("Wrong Number of Harmonics in 'GF_NSxGradHelmholtz'"); Message::Error("Wrong Number of Harmonics in 'GF_NSxGradHelmholtz'");
V->Type = SCALAR ; V->Type = SCALAR;
switch((int)Fct->Para[0]){ switch((int)Fct->Para[0]) {
case DIM_2D : case DIM_2D:
xxs = Current.x-Current.xs ; xxs = Current.x - Current.xs;
yys = Current.y-Current.ys ; yys = Current.y - Current.ys;
r = sqrt(SQU(xxs)+SQU(yys)) ; r = sqrt(SQU(xxs) + SQU(yys));
if(Current.Element->Num == NO_ELEMENT) if(Current.Element->Num == NO_ELEMENT)
Current.Element = Current.ElementSource ; Current.Element = Current.ElementSource;
ny = - Current.Element->x[1] + Current.Element->x[0] ; ny = -Current.Element->x[1] + Current.Element->x[0];
nx = Current.Element->y[1] - Current.Element->y[0] ; nx = Current.Element->y[1] - Current.Element->y[0];
n = sqrt(SQU(nx)+SQU(ny)) ; n = sqrt(SQU(nx) + SQU(ny));
nx = nx / n ; nx = nx / n;
ny = ny / n ; ny = ny / n;
if (!r) Cal_ZeroValue(V); if(!r)
Cal_ZeroValue(V);
else { else {
c1 = Fct->Para[1]/4/r ; c1 = Fct->Para[1] / 4 / r;
cr = c1 * y1(Fct->Para[1]*r); cr = c1 * y1(Fct->Para[1] * r);
ci = c1 * j1(Fct->Para[1]*r); ci = c1 * j1(Fct->Para[1] * r);
V->Val[0] = nx * xxs * cr + ny * yys * cr ; V->Val[0] = nx * xxs * cr + ny * yys * cr;
V->Val[MAX_DIM ] = nx * xxs * ci + ny * yys * ci ; V->Val[MAX_DIM] = nx * xxs * ci + ny * yys * ci;
} }
break ; break;
case DIM_3D : case DIM_3D:
xxs = Current.x-Current.xs ; xxs = Current.x - Current.xs;
yys = Current.y-Current.ys ; yys = Current.y - Current.ys;
zzs = Current.z-Current.zs ; zzs = Current.z - Current.zs;
r = sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)) ; r = sqrt(SQU(xxs) + SQU(yys) + SQU(zzs));
if (!r) Cal_ZeroValue(V); if(!r)
Cal_ZeroValue(V);
else { else {
x1x0 = Current.Element->x[1] - Current.Element->x[0] ; x1x0 = Current.Element->x[1] - Current.Element->x[0];
y1y0 = Current.Element->y[1] - Current.Element->y[0] ; y1y0 = Current.Element->y[1] - Current.Element->y[0];
z1z0 = Current.Element->z[1] - Current.Element->z[0] ; z1z0 = Current.Element->z[1] - Current.Element->z[0];
x2x0 = Current.Element->x[2] - Current.Element->x[0] ; x2x0 = Current.Element->x[2] - Current.Element->x[0];
y2y0 = Current.Element->y[2] - Current.Element->y[0] ; y2y0 = Current.Element->y[2] - Current.Element->y[0];
z2z0 = Current.Element->z[2] - Current.Element->z[0] ; z2z0 = Current.Element->z[2] - Current.Element->z[0];
nx = y1y0 * z2z0 - z1z0 * y2y0 ; nx = y1y0 * z2z0 - z1z0 * y2y0;
ny = z1z0 * x2x0 - x1x0 * z2z0 ; ny = z1z0 * x2x0 - x1x0 * z2z0;
nz = x1x0 * y2y0 - y1y0 * x2x0 ; nz = x1x0 * y2y0 - y1y0 * x2x0;
n = sqrt(SQU(nx)+SQU(ny)+SQU(nz)) ; n = sqrt(SQU(nx) + SQU(ny) + SQU(nz));
nx = nx/n ; nx = nx / n;
ny = ny/n ; ny = ny / n;
nz = nz/n ; nz = nz / n;
c1 = - ONE_OVER_FOUR_PI / CUB(r) ; c1 = -ONE_OVER_FOUR_PI / CUB(r);
c2 = ONE_OVER_FOUR_PI * Fct->Para[1] / SQU(r) ; c2 = ONE_OVER_FOUR_PI * Fct->Para[1] / SQU(r);
cr = (c1 * cos(Fct->Para[1]*r) - c2 * sin(Fct->Para[1]*r)) ; cr = (c1 * cos(Fct->Para[1] * r) - c2 * sin(Fct->Para[1] * r));
ci = (c1 * sin(Fct->Para[1]*r) + c2 * cos(Fct->Para[1]*r)) ; ci = (c1 * sin(Fct->Para[1] * r) + c2 * cos(Fct->Para[1] * r));
V->Val[0] =nx * xxs * cr + ny * yys * cr + nz * zzs * cr ; V->Val[0] = nx * xxs * cr + ny * yys * cr + nz * zzs * cr;
V->Val[MAX_DIM ] = nx* xxs * ci + ny * yys * ci + nz * zzs * ci; V->Val[MAX_DIM] = nx * xxs * ci + ny * yys * ci + nz * zzs * ci;
} }
break ; break;
default : default:
Message::Error("Bad Parameter for 'GF_NSxGradHelmholtz' (%d)", (int)Fct->Par Message::Error("Bad Parameter for 'GF_NSxGradHelmholtz' (%d)",
a[0]); (int)Fct->Para[0]);
break; break;
} }
} }
 End of changes. 47 change blocks. 
170 lines changed or deleted 174 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)