is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| Bessel.cpp (getdp-3.4.0-source.tgz) | : | Bessel.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. |
| | |
| #include "GetDPConfig.h" | | #include "GetDPConfig.h" |
| #include "Message.h" | | #include "Message.h" |
| #include "Bessel.h" | | #include "Bessel.h" |
| | |
| #if defined(HAVE_NO_FORTRAN) | | #if defined(HAVE_NO_FORTRAN) |
| | |
|
| static void zbesj_(double*, double*, double*, int*, int*, double*, | | static void zbesj_(double *, double *, double *, int *, int *, double *, |
| double*, int*, int*) | | double *, int *, int *) |
| { | | { |
| Message::Fatal("Bessel functions require Fortran compiler"); | | Message::Fatal("Bessel functions require Fortran compiler"); |
| } | | } |
| | |
|
| static void zbesk_(double*, double*, double*, int*, int*, double*, | | static void zbesk_(double *, double *, double *, int *, int *, double *, |
| double*, int*, int*) | | double *, int *, int *) |
| { | | { |
| Message::Fatal("Bessel functions require Fortran compiler"); | | Message::Fatal("Bessel functions require Fortran compiler"); |
| } | | } |
| | |
|
| static void zbesy_(double*, double*, double*, int*, int*, double*, | | static void zbesy_(double *, double *, double *, int *, int *, double *, |
| double*, int*, double*, double*, int*) | | double *, int *, double *, double *, int *) |
| { | | { |
| Message::Fatal("Bessel functions require Fortran compiler"); | | Message::Fatal("Bessel functions require Fortran compiler"); |
| } | | } |
| | |
|
| static void zbesh_(double*, double*, double*, int*, int*, int*, | | static void zbesh_(double *, double *, double *, int *, int *, int *, double *, |
| double*, double*, int*, int*) | | double *, int *, int *) |
| { | | { |
| Message::Fatal("Bessel functions require Fortran compiler"); | | Message::Fatal("Bessel functions require Fortran compiler"); |
| } | | } |
| | |
| #else | | #else |
| | |
| #if defined(HAVE_UNDERSCORE) | | #if defined(HAVE_UNDERSCORE) |
| #define zbesj_ zbesj | | #define zbesj_ zbesj |
| define zbesk_ zbesk | | define zbesk_ zbesk |
| #define zbesy_ zbesy | | #define zbesy_ zbesy |
| #define zbesh_ zbesh | | #define zbesh_ zbesh |
| #endif | | #endif |
| | |
|
| extern "C" { | | extern "C" |
| void zbesj_(double*, double*, double*, int*, int*, double*, | | { |
| double*, int*, int*); | | void zbesj_(double *, double *, double *, int *, int *, double *, double *, |
| void zbesk_(double*, double*, double*, int*, int*, double*, | | int *, int *); |
| double*, int*, int*); | | void zbesk_(double *, double *, double *, int *, int *, double *, double *, |
| | | int *, int *); |
| void zbesy_(double*, double*, double*, int*, int*, double*, | | |
| double*, int*, double*, | | void zbesy_(double *, double *, double *, int *, int *, double *, double *, |
| double*, int*); | | int *, double *, double *, int *); |
| void zbesh_(double*, double*, double*, int*, int*, int*, double*, | | void zbesh_(double *, double *, double *, int *, int *, int *, double *, |
| double*, int*, int*); | | double *, int *, int *); |
| } | | } |
| | |
| #endif | | #endif |
| | |
| static int BesselError(int ierr, const char *str) | | static int BesselError(int ierr, const char *str) |
| { | | { |
|
| static int warn=0; | | static int warn = 0; |
| | |
|
| switch(ierr){ | | switch(ierr) { |
| case 0 : | | case 0: return 0; |
| return 0; | | case 1: Message::Error("Input error in %s", str); return BESSEL_ERROR_INPUT; |
| case 1 : | | case 2: return BESSEL_OVERFLOW; |
| Message::Error("Input error in %s", str); | | case 3: |
| return BESSEL_ERROR_INPUT; | | if(!warn) { |
| case 2 : | | Message::Info( |
| return BESSEL_OVERFLOW; | | "Half machine accuracy lost in %s (large argument or order)", str); |
| case 3 : | | |
| if(!warn){ | | |
| Message::Info("Half machine accuracy lost in %s (large argument or order)" | | |
| , str); | | |
| warn = 1; | | warn = 1; |
| } | | } |
| return BESSEL_HALF_ACCURACY; | | return BESSEL_HALF_ACCURACY; |
|
| case 4 : | | case 4: |
| Message::Error("Complete loss of significance in %s (argument or order too l | | Message::Error( |
| arge)", str); | | "Complete loss of significance in %s (argument or order too large)", str); |
| return BESSEL_NO_ACCURACY; | | return BESSEL_NO_ACCURACY; |
|
| case 5 : | | case 5: |
| Message::Error("Failed to converge in %s", str); | | Message::Error("Failed to converge in %s", str); |
| return BESSEL_NO_CONVERGENCE; | | return BESSEL_NO_CONVERGENCE; |
| default: | | default: |
| Message::Info("Unknown Bessel status in %s (%d)", str, ierr); | | Message::Info("Unknown Bessel status in %s (%d)", str, ierr); |
| return ierr; | | return ierr; |
| } | | } |
| } | | } |
| | |
| // First kind Bessel functions | | // First kind Bessel functions |
| | |
| int BesselJn(double n, int num, double x, double *val) | | int BesselJn(double n, int num, double x, double *val) |
| { | | { |
| int nz = 0, ierr = 0, kode = 1; | | int nz = 0, ierr = 0, kode = 1; |
| double xi = 0.0; | | double xi = 0.0; |
|
| double* ji = new double[num]; | | double *ji = new double[num]; |
| | |
|
| zbesj_(&x, &xi, &n, &kode, &num, val, ji, &nz, &ierr) ; | | zbesj_(&x, &xi, &n, &kode, &num, val, ji, &nz, &ierr); |
| | |
| delete[] ji; | | delete[] ji; |
| | |
| return BesselError(ierr, "BesselJn"); | | return BesselError(ierr, "BesselJn"); |
| } | | } |
| | |
|
| int BesselJnComplex(double n, int num, double xr, double xi, double *valr, doubl | | int BesselJnComplex(double n, int num, double xr, double xi, double *valr, |
| e *vali) | | double *vali) |
| { | | { |
| int nz = 0, ierr = 0, kode = 1; | | int nz = 0, ierr = 0, kode = 1; |
| | |
|
| zbesj_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr) ; | | zbesj_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr); |
| | |
| return BesselError(ierr, "BesselJnComplex"); | | return BesselError(ierr, "BesselJnComplex"); |
| } | | } |
| | |
|
| int BesselKnComplex(double n, int num, double xr, double xi, double *valr, doubl | | int BesselKnComplex(double n, int num, double xr, double xi, double *valr, |
| e *vali) | | double *vali) |
| { | | { |
| int nz = 0, ierr = 0, kode = 1; | | int nz = 0, ierr = 0, kode = 1; |
| | |
|
| zbesk_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr) ; | | zbesk_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr); |
| | |
| return BesselError(ierr, "BesselKnComplex"); | | return BesselError(ierr, "BesselKnComplex"); |
| } | | } |
| | |
| int BesselSphericalJn(double n, int num, double x, double *val) | | int BesselSphericalJn(double n, int num, double x, double *val) |
| { | | { |
|
| int ierr = BesselJn(n+0.5, num, x, val); | | int ierr = BesselJn(n + 0.5, num, x, val); |
| double coef = sqrt(0.5*M_PI/x); | | double coef = sqrt(0.5 * M_PI / x); |
| for(int i = 0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] *= coef; } |
| val[i] *= coef; | | |
| } | | |
| return BesselError(ierr, "BesselSphericalJn"); | | return BesselError(ierr, "BesselSphericalJn"); |
| } | | } |
| | |
| int BesselAltSphericalJn(double n, int num, double x, double *val) | | int BesselAltSphericalJn(double n, int num, double x, double *val) |
| { | | { |
|
| int ierr = BesselJn(n+0.5, num, x, val); | | int ierr = BesselJn(n + 0.5, num, x, val); |
| double coef = sqrt(0.5*M_PI*x); | | double coef = sqrt(0.5 * M_PI * x); |
| for(int i = 0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] *= coef; } |
| val[i] *= coef; | | |
| } | | |
| return BesselError(ierr, "BesselAltSphericalJn"); | | return BesselError(ierr, "BesselAltSphericalJn"); |
| } | | } |
| | |
| // Second kind Bessel functions | | // Second kind Bessel functions |
| | |
| int BesselYn(double n, int num, double x, double *val) | | int BesselYn(double n, int num, double x, double *val) |
| { | | { |
| int nz = 0, ierr = 0, kode = 1; | | int nz = 0, ierr = 0, kode = 1; |
| double xi = 0.0; | | double xi = 0.0; |
|
| double* yi = new double[num]; | | double *yi = new double[num]; |
| double* auxyr = new double[num]; | | double *auxyr = new double[num]; |
| double* auxyi = new double[num]; | | double *auxyi = new double[num]; |
| | |
| zbesy_(&x, &xi, &n, &kode, &num, val, yi, &nz, auxyr, auxyi, &ierr); | | zbesy_(&x, &xi, &n, &kode, &num, val, yi, &nz, auxyr, auxyi, &ierr); |
| | |
| delete[] yi; | | delete[] yi; |
| delete[] auxyr; | | delete[] auxyr; |
| delete[] auxyi; | | delete[] auxyi; |
| | |
| return BesselError(ierr, "BesselYn"); | | return BesselError(ierr, "BesselYn"); |
| } | | } |
| | |
| int BesselSphericalYn(double n, int num, double x, double *val) | | int BesselSphericalYn(double n, int num, double x, double *val) |
| { | | { |
|
| int ierr = BesselYn(n+0.5, num, x, val); | | int ierr = BesselYn(n + 0.5, num, x, val); |
| double coef = sqrt(0.5*M_PI/x); | | double coef = sqrt(0.5 * M_PI / x); |
| for(int i = 0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] *= coef; } |
| val[i] *= coef; | | |
| } | | |
| return BesselError(ierr, "BesselSphericalYn"); | | return BesselError(ierr, "BesselSphericalYn"); |
| } | | } |
| | |
| int BesselAltSphericalYn(double n, int num, double x, double *val) | | int BesselAltSphericalYn(double n, int num, double x, double *val) |
| { | | { |
|
| int ierr = BesselYn(n+0.5, num, x, val); | | int ierr = BesselYn(n + 0.5, num, x, val); |
| double coef = sqrt(0.5*M_PI*x); | | double coef = sqrt(0.5 * M_PI * x); |
| for(int i = 0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] *= coef; } |
| val[i] *= coef; | | |
| } | | |
| return BesselError(ierr, "BesselAltSphericalYn"); | | return BesselError(ierr, "BesselAltSphericalYn"); |
| } | | } |
| | |
| // Hankel functions (type = 1 or 2) | | // Hankel functions (type = 1 or 2) |
| | |
| int BesselHn(int type, double n, int num, double x, std::complex<double> *val) | | int BesselHn(int type, double n, int num, double x, std::complex<double> *val) |
| { | | { |
| int nz = 0, ierr = 0, kode = 1; | | int nz = 0, ierr = 0, kode = 1; |
|
| double* hr = new double[num]; | | double *hr = new double[num]; |
| double* hi = new double[num]; | | double *hi = new double[num]; |
| double xi = 0.0; | | double xi = 0.0; |
| | |
| zbesh_(&x, &xi, &n, &kode, &type, &num, hr, hi, &nz, &ierr); | | zbesh_(&x, &xi, &n, &kode, &type, &num, hr, hi, &nz, &ierr); |
| | |
|
| for(int i=0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] = std::complex<double>(hr[i], hi[i]); } |
| val[i] = std::complex<double>(hr[i], hi[i]); | | |
| } | | |
| | |
| delete[] hr; | | delete[] hr; |
| delete[] hi; | | delete[] hi; |
| | |
| return BesselError(ierr, "BesselHn"); | | return BesselError(ierr, "BesselHn"); |
| } | | } |
| | |
|
| int BesselSphericalHn(int type, double n, int num, double x, std::complex<double | | int BesselSphericalHn(int type, double n, int num, double x, |
| > *val) | | std::complex<double> *val) |
| { | | { |
|
| int ierr = BesselHn(type, n+0.5, num, x, val); | | int ierr = BesselHn(type, n + 0.5, num, x, val); |
| double coef = sqrt(0.5*M_PI/x); | | double coef = sqrt(0.5 * M_PI / x); |
| for(int i = 0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] *= coef; } |
| val[i] *= coef; | | |
| } | | |
| return BesselError(ierr, "BesselSphericalHn"); | | return BesselError(ierr, "BesselSphericalHn"); |
| } | | } |
| | |
|
| int BesselAltSphericalHn(int type, double n, int num, double x, std::complex<dou | | int BesselAltSphericalHn(int type, double n, int num, double x, |
| ble> *val) | | std::complex<double> *val) |
| { | | { |
|
| int ierr = BesselHn(type, n+0.5, num, x, val); | | int ierr = BesselHn(type, n + 0.5, num, x, val); |
| double coef = sqrt(0.5*M_PI*x); | | double coef = sqrt(0.5 * M_PI * x); |
| for(int i = 0; i < num; i++){ | | for(int i = 0; i < num; i++) { val[i] *= coef; } |
| val[i] *= coef; | | |
| } | | |
| return BesselError(ierr, "BesselAltSphericalHn"); | | return BesselError(ierr, "BesselAltSphericalHn"); |
| } | | } |
| | |
| // Utilities for backward compatibility | | // Utilities for backward compatibility |
| | |
| double Spherical_j_n(int n, double x) | | double Spherical_j_n(int n, double x) |
| { | | { |
| double res; | | double res; |
| BesselSphericalJn(n, 1, x, &res); | | BesselSphericalJn(n, 1, x, &res); |
| return res; | | return res; |
| | | |
| End of changes. 27 change blocks. |
|---|
| 87 lines changed or deleted | | 69 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 