diff --git a/contrib/FourierModel/Interpolator1D.cpp b/contrib/FourierModel/Interpolator1D.cpp
deleted file mode 100755
index 9d51584f438ab64ae0abcab27ea2442687166229..0000000000000000000000000000000000000000
--- a/contrib/FourierModel/Interpolator1D.cpp
+++ /dev/null
@@ -1,602 +0,0 @@
-#include <vector>
-#include <math.h>
-#include <complex>
-#include "Interpolator1D.h"
-
-FftPolyInterpolator1D::FftPolyInterpolator1D(const std::vector<double> &u,
- const std::vector< std::complex<double> > &data,
- int refineFactor, int polyOrder, int derivative,
- bool storeCoefs)
- : _refineFactor(refineFactor), _polyOrder(polyOrder), _derivative(derivative),
- _storeCoefs(storeCoefs), _uFine(0), _fineData(0), _fineDeriv(0), _fineDeriv2(0)
-{
- _uIntervals = u.size() - 1;
- _uFineIntervals = _refineFactor * _uIntervals;
-
- // min/max of input mesh
- if(u[_uIntervals] > u[0]) {
- _uMin = u[0];
- _uMax = u[_uIntervals];
- }
- else {
- _uMin = u[_uIntervals];
- _uMax = u[0];
- }
-
- // signed length of interval
- _uLength = u[_uIntervals] - u[0];
-
- // Sanity checks
- if(std::abs(data[0] - data[_uIntervals]) > 1.e-10){
- Msg::Warning("Data is not periodic: f(%g)=(%.16g,%.16g), f(%g)=(%.16g,%.16g)",
- u[0], data[0].real(), data[0].imag(),
- u[_uIntervals], data[_uIntervals].real(), data[_uIntervals].imag());
- }
- if(fabs(u[1] - u[0]) - fabs(u[2] - u[1]) > 1.e-10){
- Msg::Fatal("Input grid is not equispaced");
- }
- if(_uFineIntervals < _polyOrder){
- Msg::Fatal("Number of intervals is smaller than polynomial order (%d < %d)",
- _uFineIntervals, _polyOrder);
- }
-
- // Compute the fine grid spacing (assuming the u grid is equally spaced)
- _hUFine = (u[1] - u[0]) / _refineFactor;
-
- // Compute the points on the fine grid
- _uFine = new double[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++)
- _uFine[i] = u[0] + i * _hUFine;
-
- // Initialize _fineData to zero
- _fineData = new std::complex<double>[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++)
- _fineData[i] = 0.;
-
- if(_storeCoefs){
- // number of local interpolating polynomials
- _nLocPol = 1;
- int flag = 0;
- while ((flag + _polyOrder) < (_uFineIntervals + 1)) {
- _nLocPol++;
- flag += _polyOrder;
- }
- if ((_uFineIntervals + 1 - flag) == 1) _nLocPol--;
-
- // local interpolating polynomial coefficients
- _locpol = new std::complex<double>[_nLocPol * (_polyOrder + 1)];
- for(int i = 0; i < _nLocPol * (_polyOrder + 1); i++)
- _locpol[i] = 0.;
- }
-
- std::complex<double> *forwardData = new std::complex<double>[_uIntervals+1];
-
- if (_refineFactor == 1){ // No Fourier interpolation
- // Copy data into _fineData
- for(int i = 0; i < _uFineIntervals + 1; i++)
- _fineData[i] = data[i];
- }
- else {
- // Compute the forward FFT of the one-dimensional data (only uses
- // data in i=0,...,uIntervals-1 and scales by 1/uIntervals)
- for(int i = 0; i < _uIntervals + 1; i++)
- forwardData[i] = data[i];
- _ForwardFft(_uIntervals + 1, forwardData);
-
- // Copy the Fourier coefficients into _fineData
- int halfIndexU = _uIntervals / 2;
- for(int i = 0; i < _uIntervals; i++){
- int iFine = i;
- if(i > halfIndexU)
- iFine = _uFineIntervals + (i - _uIntervals);
- _fineData[iFine] = forwardData[i];
- }
-
- // Compute inverse FFT on _fineData to interpolate the data on the
- // fine grid
- _BackwardFft(_uFineIntervals + 1, _fineData);
- }
-
- // Check if we need to interpolate the derivative(s)
- if(_derivative){
- if(_refineFactor == 1){
- // No Fourier interpolation was performed, so we need to do it here...
- for(int i = 0; i < _uIntervals + 1; i++)
- forwardData[i] = data[i];
- _ForwardFft(_uIntervals + 1, forwardData);
- }
-
- // Initialize _fineDeriv and _fineDeriv2 to zero
- if(_derivative & 1){
- _fineDeriv = new std::complex<double>[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++)
- _fineDeriv[i] = 0.;
- if(_storeCoefs){
- // local interpolating polynomial coefficients
- _locpolDeriv = new std::complex<double>[_nLocPol*(_polyOrder+1)];
- for(int i = 0; i < _nLocPol*(_polyOrder+1); i++)
- _locpolDeriv[i] = 0.;
- }
- }
- if(_derivative & 2){
- _fineDeriv2 = new std::complex<double>[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++)
- _fineDeriv2[i] = 0.;
- if(_storeCoefs){
- // local interpolating polynomial coefficients
- _locpolDeriv2 = new std::complex<double>[_nLocPol*(_polyOrder+1)];
- for(int i = 0; i < _nLocPol*(_polyOrder+1); i++)
- _locpolDeriv2[i] = 0.;
- }
- }
-
- // Copy the Fourier coefficients into _fineDeriv and _fineDeriv2
- std::complex<double> I(0., 1.);
- int halfIndexU = _uIntervals / 2;
- for(int i = 0; i < _uIntervals; i++){
- int iFine = i;
- double k = i;
- if(i > halfIndexU){
- k = i - _uIntervals;
- iFine = _uFineIntervals + (i - _uIntervals);
- }
- // multiply by (2*pi*i*k) to get the derivative
- if(_derivative & 1)
- _fineDeriv[iFine] = (2. * M_PI * k * I) * forwardData[i] / _uLength;
- // multiply by (2*pi*i*k)^2 to get the second derivative
- if(_derivative & 2)
- _fineDeriv2[iFine] =
- -(4. * M_PI * M_PI * k * k) * forwardData[i] / (_uLength * _uLength);
- }
-
- // Perform an inverse FFT on _fineDeriv to interpolate the
- // derivative to the fine grid
- if(_derivative & 1)
- _BackwardFft(_uFineIntervals + 1, _fineDeriv);
- if(_derivative & 2)
- _BackwardFft(_uFineIntervals + 1, _fineDeriv2);
- }
-
- delete [] forwardData;
-
- if(_storeCoefs){
- _interPolyCoeff(_locpol, _fineData, _uFineIntervals, _hUFine, _polyOrder);
- if(_derivative & 1)
- _interPolyCoeff(_locpolDeriv, _fineDeriv, _uFineIntervals, _hUFine,
- _polyOrder);
- if(_derivative & 2)
- _interPolyCoeff(_locpolDeriv2, _fineDeriv2, _uFineIntervals, _hUFine,
- _polyOrder);
- }
-
- // Initialize temp interpolation variables
- _tmpInterp = new double[_polyOrder + 1];
- _tmpPnts = new double[_polyOrder + 1];
- _tmpValsReal = new double[_polyOrder + 1];
- _tmpValsImag = new double[_polyOrder + 1];
-}
-
-FftPolyInterpolator1D::~FftPolyInterpolator1D()
-{
- delete[] _uFine;
- delete[] _fineData;
- if(_storeCoefs)
- delete[] _locpol;
- if(_fineDeriv) {
- delete[] _fineDeriv;
- if(_storeCoefs)
- delete[] _locpolDeriv;
- }
- if(_fineDeriv2) {
- delete[] _fineDeriv2;
- if(_storeCoefs)
- delete[] _locpolDeriv2;
- }
-
- delete[] _tmpInterp;
- delete[] _tmpPnts;
- delete[] _tmpValsReal;
- delete[] _tmpValsImag;
-}
-
-int FftPolyInterpolator1D::_forwardSize = 0;
-int FftPolyInterpolator1D::_backwardSize = 0;
-fftw_plan FftPolyInterpolator1D::_forwardPlan;
-fftw_plan FftPolyInterpolator1D::_backwardPlan;
-fftw_complex *FftPolyInterpolator1D::_forwardData = 0;
-fftw_complex *FftPolyInterpolator1D::_backwardData = 0;
-
-void FftPolyInterpolator1D::_SetForwardPlan(int n)
-{
- if(n != _forwardSize){
- if(_forwardSize){
- fftw_destroy_plan(_forwardPlan);
- fftw_free(_forwardData);
- }
- _forwardSize = n;
- _forwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _forwardSize);
- _forwardPlan = fftw_plan_dft_1d(_forwardSize, _forwardData, _forwardData,
- FFTW_FORWARD, FFTW_ESTIMATE);
- }
-}
-
-void FftPolyInterpolator1D::_SetBackwardPlan(int n)
-{
- if(n != _backwardSize){
- if(_backwardSize){
- fftw_destroy_plan(_backwardPlan);
- fftw_free(_backwardData);
- }
- _backwardSize = n;
- _backwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _backwardSize);
- _backwardPlan = fftw_plan_dft_1d(_backwardSize, _backwardData, _backwardData,
- FFTW_BACKWARD, FFTW_ESTIMATE);
- }
-}
-
-void FftPolyInterpolator1D::_ForwardFft(int n, std::complex<double> *fftData)
-{
- // Initialize fftw plan and array (ignoring the last element of
- // fftData, which should just be the periodic extension)
- _SetForwardPlan(n - 1);
- for(int i = 0; i < n - 1; i++){
- _forwardData[i][0] = fftData[i].real();
- _forwardData[i][1] = fftData[i].imag();
- }
-
- // Perform forward FFT
- fftw_execute(_forwardPlan);
-
- // Copy data back into fftData and scale by 1/(n - 1)
- double s = 1. / (double)(n - 1);
- for(int i = 0; i < n - 1; i++)
- fftData[i] = s * std::complex<double>(_forwardData[i][0], _forwardData[i][1]);
-}
-
-void FftPolyInterpolator1D::_BackwardFft(int n, std::complex<double> *fftData)
-{
- // Initialize fftw plan and array (ignoring last element of fftData)
- _SetBackwardPlan(n - 1);
- for(int i = 0; i < n - 1; i++){
- _backwardData[i][0] = fftData[i].real();
- _backwardData[i][1] = fftData[i].imag();
- }
-
- // Perform backward FFT
- fftw_execute(_backwardPlan);
-
- // Copy data back into fftData
- for(int i = 0; i < n - 1; i++)
- fftData[i] = std::complex<double>(_backwardData[i][0], _backwardData[i][1]);
-
- // Fill in last element with copy of first element
- fftData[n - 1] = fftData[0];
-}
-
-double FftPolyInterpolator1D::_PolyInterp(double *pnts, double *vals, double t)
-{
- // Neville's Algorithm from Stoer and Bulirsch, Section 2.1.2
- for(int i = 0; i <= _polyOrder; i++){
- _tmpInterp[i] = vals[i];
- for(int j = i-1; j >= 0; j--)
- _tmpInterp[j] = _tmpInterp[j+1] +
- (_tmpInterp[j+1] - _tmpInterp[j]) * (t - pnts[i]) / (pnts[i] - pnts[j]);
- }
- return _tmpInterp[0];
-}
-
-std::complex<double> FftPolyInterpolator1D::_Interpolate(double u, int derivative)
-{
- if((derivative == 1 && !(_derivative & 1)) ||
- (derivative == 2 && !(_derivative & 2)) ||
- derivative < 0 || derivative > 2){
- Msg::Error("Derivative data not available: check constructor call");
- return 0.;
- }
-
- if(u == _uMax)
- u = _uMin;
-
- double epsilon = 1e-12;
- if(u < _uMin - epsilon || u > _uMax + epsilon){
- Msg::Error("Trying to interpolate outside interval: u=%.16g not in [%g,%g]",
- u, _uMin, _uMax);
- return 0.;
- }
-
- // Choose uIndexStart so that u is centered in the polynomial
- // interpolation grid
- int uIndexStart = (int)floor((u - _uFine[0]) / _hUFine) - _polyOrder / 2;
- if(uIndexStart < 0)
- uIndexStart = 0;
- else if((uIndexStart + _polyOrder) > _uFineIntervals)
- uIndexStart = _uFineIntervals - _polyOrder;
-
- // Interpolate to find value at u
- for(int i = 0; i <= _polyOrder; i++){
- _tmpPnts[i] = _uFine[uIndexStart + i];
- if(derivative == 0){
- _tmpValsReal[i] = _fineData[uIndexStart + i].real();
- _tmpValsImag[i] = _fineData[uIndexStart + i].imag();
- }
- else if(derivative == 1){
- _tmpValsReal[i] = _fineDeriv[uIndexStart + i].real();
- _tmpValsImag[i] = _fineDeriv[uIndexStart + i].imag();
- }
- else{
- _tmpValsReal[i] = _fineDeriv2[uIndexStart + i].real();
- _tmpValsImag[i] = _fineDeriv2[uIndexStart + i].imag();
- }
- }
-
- return std::complex<double>(_PolyInterp(_tmpPnts, _tmpValsReal, u),
- _PolyInterp(_tmpPnts, _tmpValsImag, u));
-}
-
-void polyint(double *xa, double *ya, int n, double x, double &y)
-{
-#define ind(ii) ((ii) - 1)
- int i, m, ns = 1;
- double dy, den, dif, dift, ho, hp, w;
-
- double *c = new double[n];
- double *d = new double[n];
-
- dif = fabs(x - xa[ind(1)]);
- for (i = 1; i <= n; i++) {
- if ((dift = fabs(x - xa[ind(i)])) < dif) {
- ns = i;
- dif = dift;
- }
- c[ind(i)] = ya[ind(i)];
- d[ind(i)] = ya[ind(i)];
- }
- y = ya[ind(ns--)];
- for (m = 1; m < n; m++) {
- for (i = 1; i <= n - m; i++) {
- ho = xa[ind(i)] - x;
- hp = xa[ind(i + m)] - x;
- w = c[ind(i + 1)] - d[ind(i)];
- if ((den = ho - hp) == 0) {
- Msg::Error("Error in POLYINT");
- return;
- }
- den = w / den;
- d[ind(i)] = hp * den;
- c[ind(i)] = ho * den;
- }
- if (2 * ns < (n - m))
- dy = c[ind(ns + 1)];
- else {
- dy = d[ind(ns)];
- ns--;
- }
- y += dy;
- }
- delete[] c;
- delete[] d;
-}
-
-void polcof(double *x, double *y, int n, double *cof)
-{
- int k;
- double xmin;
-
- double *xtemp = new double[n + 1];
- double *ytemp = new double[n + 1];
-
- for (int j = 0; j <= n; j++) {
- xtemp[j] = x[j];
- ytemp[j] = y[j];
- }
- for (int j = 0; j <= n; j++) {
- polyint(xtemp, ytemp, n + 1 - j, 0.0, cof[j]);
- xmin = 1.0e38;
- k = -1;
- for (int i = 0; i <= n - j; i++) {
- if (fabs(xtemp[i]) < xmin) {
- xmin = fabs(xtemp[i]);
- k = i;
- }
- if (xtemp[i] != 0) ytemp[i] = (ytemp[i] - cof[j]) / xtemp[i];
- }
- for (int i = k + 1; i<= n - j; i++) {
- xtemp[i - 1] = xtemp[i];
- ytemp[i - 1] = ytemp[i];
- }
- }
- delete[] xtemp;
- delete[] ytemp;
-}
-
-void getcoeff(double *f,double *x,int n)
-{
- double *ya = new double[n];
- double *cof = new double[n];
-
- for (int i = 0; i < n; i++) {
- ya[i] = f[i];
- }
- polcof(x, ya, n - 1, cof);
- for (int i = 0; i < n; i++)
- f[i] = cof[i];
-
- delete[] ya;
- delete[] cof;
-}
-
-void FftPolyInterpolator1D::_interPolyCoeff(std::complex<double> *locpol,
- std::complex<double> *u,
- int rn, double fineh, int npoly)
-{
- double *x = new double[npoly + 1];
- double *freal = new double[npoly + 1];
- double *fimag = new double[npoly + 1];
- int rem;
-
- int l = 0;
- int flag = 0;
- while ((flag + npoly) < rn) {
- for (int i = 0; i <= npoly; i++)
- x[i] = (flag + i) * fineh;
- for (int i = 0; i <= npoly; i++) {
- freal[i] = real(u[flag + i]);
- fimag[i] = imag(u[flag + i]);
- }
- getcoeff(freal, x, npoly + 1);
- getcoeff(fimag, x, npoly + 1);
- for (int i = 0; i < npoly + 1; i++)
- locpol[(npoly + 1) * l + i] = std::complex<double>(freal[i], fimag[i]);
- l++;
- flag += npoly;
- }
-
- rem = rn - flag;
- if (rem > 1) {
- for (int i = 0; i <= npoly; i++)
- x[i] = (flag + i - (npoly + 1) + rem) * fineh;
- for (int i=0;i<=npoly;i++) {
- freal[i] = real(u[flag + i - (npoly + 1) + rem]);
- fimag[i] = imag(u[flag + i - (npoly + 1) + rem]);
- }
- getcoeff(freal, x, npoly + 1);
- getcoeff(fimag, x, npoly + 1);
- for (int i = 0; i < npoly + 1; i++)
- locpol[(npoly + 1) * l + i] = std::complex<double>(freal[i], fimag[i]);
- }
-
- delete[] x;
- delete[] freal;
- delete[] fimag;
-}
-
-std::complex<double> FftPolyInterpolator1D::_InterpolateFromCoeffs(double u,
- int derivative)
-{
- if((derivative == 1 && !(_derivative & 1)) ||
- (derivative == 2 && !(_derivative & 2)) ||
- derivative < 0 || derivative > 2){
- Msg::Error("Derivative data not available: check constructor call");
- return 0.;
- }
-
- if(u == _uMax)
- u = _uMin;
-
- double epsilon = 1e-12;
- if(u < _uMin - epsilon || u > _uMax + epsilon){
- Msg::Error("Trying to interpolate outside interval: u=%.16g not in [%g,%g]",
- u, _uMin, _uMax);
- return 0.;
- }
-
- std::complex<double> out(0.,0.);
-
- // find the polynomial to be used for interpolating at u
- int l = (int)floor((u - _uFine[0])/(_polyOrder * _hUFine));
-
- if (derivative == 0) {
- out = u * _locpol[(_polyOrder + 1) * l + _polyOrder];
- for (int r = _polyOrder - 1; r > 0; r--)
- out = u * (out + _locpol[(_polyOrder + 1) * l + r]);
- out = out + _locpol[(_polyOrder + 1) * l + 0];
- }
- else if (derivative == 1) {
- out = u * _locpolDeriv[(_polyOrder + 1) * l + _polyOrder];
- for (int r = _polyOrder - 1; r > 0 ; r--)
- out = u * (out + _locpolDeriv[(_polyOrder + 1) * l + r]);
- out = out + _locpolDeriv[(_polyOrder + 1) * l + 0];
- }
- else if (derivative == 2) {
- out = u * _locpolDeriv2[(_polyOrder + 1) * l + _polyOrder];
- for (int r = _polyOrder - 1; r > 0; r--)
- out = u * (out + _locpolDeriv2[(_polyOrder + 1) * l + r]);
- out = out + _locpolDeriv2[(_polyOrder + 1) * l + 0];
- }
-
- return out;
-}
-
-std::complex<double> FftPolyInterpolator1D::F(double u)
-{
- if(_storeCoefs)
- return _InterpolateFromCoeffs(u, 0);
- else
- return _Interpolate(u, 0);
-}
-
-std::complex<double> FftPolyInterpolator1D::Dfdu(double u)
-{
- if(_storeCoefs)
- return _InterpolateFromCoeffs(u, 1);
- else
- return _Interpolate(u, 1);
-}
-
-std::complex<double> FftPolyInterpolator1D::Dfdfdudu(double u)
-{
- if(_storeCoefs)
- return _InterpolateFromCoeffs(u, 2);
- else
- return _Interpolate(u, 2);
-}
-
-std::complex<double> FftPolyInterpolator1D::Integrate()
-{
- std::complex<double> value(0., 0.);
- for(int i = 0; i < _uFineIntervals; i++){
- value += (0.5 * _fineData[i] + 0.5 * _fineData[i + 1]) *
- fabs(_uFine[i + 1] - _uFine[i]);
- }
- return value;
-}
-
-FftSplineInterpolator1D::FftSplineInterpolator1D(const std::vector<double> &u,
- const std::vector< std::complex<double> > &data,
- int refineFactor)
- : FftPolyInterpolator1D(u, data, refineFactor, 3, 2)
-{
-}
-
-std::complex<double> FftSplineInterpolator1D::_SplineInterp(double *pnts,
- std::complex<double> *vals,
- std::complex<double> *deriv,
- double t)
-{
- double h = pnts[1] - pnts[0];
- double a = (pnts[1] - t) / h;
- double b = (t - pnts[0]) / h;
- double a1 = a * a * a - a;
- double b1 = b * b * b - b;
- double h1 = (h * h) / 6.;
- return a * vals[0] + b * vals[1] + (a1 * deriv[0] + b1 * deriv[1]) * h1;
-}
-
-std::complex<double> FftSplineInterpolator1D::F(double u)
-{
- double pnts[2];
- std::complex<double> vals[2], deriv[2];
-
- double epsilon = 1e-12;
- if(u < _uMin - epsilon || u > _uMax + epsilon){
- Msg::Error("Trying to interpolate outside interval: u=%.16g not in [%g,%g]",
- u, _uMin, _uMax);
- return 0.;
- }
-
- // Get indices of grid points surrounding u
- int uIndexStart = (int)floor((u - _uFine[0]) / _hUFine);
- if(uIndexStart < 0)
- uIndexStart = 0;
- else if(uIndexStart > _uFineIntervals - 1)
- uIndexStart = _uFineIntervals - 1;
-
- pnts[0] = _uFine[uIndexStart];
- pnts[1] = _uFine[uIndexStart + 1];
- vals[0] = _fineData[uIndexStart];
- vals[1] = _fineData[uIndexStart + 1];
- deriv[0] = _fineDeriv2[uIndexStart];
- deriv[1] = _fineDeriv2[uIndexStart + 1];
-
- return _SplineInterp(pnts, vals, deriv, u);
-}
diff --git a/contrib/FourierModel/Interpolator1D.h b/contrib/FourierModel/Interpolator1D.h
deleted file mode 100755
index 0db89ae1c4224aef0ae9160cdd7ba07a6c0610c9..0000000000000000000000000000000000000000
--- a/contrib/FourierModel/Interpolator1D.h
+++ /dev/null
@@ -1,104 +0,0 @@
-#ifndef _INTERPOLATOR_1D_H_
-#define _INTERPOLATOR_1D_H_
-
-#include <vector>
-#include <complex>
-#include <fftw3.h>
-#include "Message.h"
-
-// The base class for the 1D interpolators
-class Interpolator1D {
- public:
- Interpolator1D(){}
- virtual ~Interpolator1D(){}
- // interpolates the data at u
- virtual std::complex<double> F(double u) = 0;
- // interpolates the first derivative of the data at u
- virtual std::complex<double> Dfdu(double u) = 0;
- // interpolates the second derivative of the data at u
- virtual std::complex<double> Dfdfdudu(double u) = 0;
-};
-
-// FFT + polynomial interpolation on refined grid (assumes that the
-// input grid is equally spaced and the input data is periodic)
-class FftPolyInterpolator1D : public Interpolator1D {
- private:
- // refinement factor (e.g. 16; if equal to 1, we just perform
- // standard piecewise polynomial interpolation, without any FFTs)
- int _refineFactor;
- // order of the interpolating polynomial
- int _polyOrder;
- // signed length of the interval
- double _uLength;
- // temporary interpolation variables
- double *_tmpInterp, *_tmpPnts, *_tmpValsReal, *_tmpValsImag;
- // interpolation polynomial coefficients computation routine
- void _interPolyCoeff(std::complex<double> *locpol, std::complex<double> *u,
- int rn,double fineh,int npoly);
- // 1D polynomial interpolation routine
- double _PolyInterp(double *pnts, double *vals, double t);
- // interpolation wrappers
- std::complex<double> _Interpolate(double u, int derivative=0);
- std::complex<double> _InterpolateFromCoeffs(double u, int derivative=0);
- // persistent fftw data (used to avoid recomputing the FFTW plans
- // and reallocating memory every time)
- static int _forwardSize, _backwardSize;
- static fftw_plan _forwardPlan, _backwardPlan;
- static fftw_complex *_forwardData, *_backwardData;
- void _SetForwardPlan(int n);
- void _SetBackwardPlan(int n);
- protected:
- // number of intervals in the original and in the refined mesh
- int _uIntervals, _uFineIntervals;
- // fine mesh spacing
- double _hUFine;
- // min/max of input mesh
- double _uMin, _uMax;
- // bitfield telling if we also interpolate the derivative(s)
- int _derivative;
- // flag set if we should store the polynomial coefficients
- bool _storeCoefs;
- // grid points of the refined mesh
- double *_uFine;
- // data (and its first 2 derivatives) on the refined regular grid
- std::complex<double> *_fineData, *_fineDeriv, *_fineDeriv2;
- // number of local interpolating polynomials
- int _nLocPol;
- // polynomial coefficients
- std::complex<double> *_locpol, *_locpolDeriv, *_locpolDeriv2;
- // This routine computes the forward FFT (ignoring the last element
- // of the data)
- void _ForwardFft(int n, std::complex<double> *fftData);
- // This routine computes the inverse FFT (ignoring the last element
- // in the array), returns values in entire array (by using the
- // periodic extension)
- void _BackwardFft(int n, std::complex<double> *fftData);
- public:
- FftPolyInterpolator1D(const std::vector<double> &u,
- const std::vector< std::complex<double> > &data,
- int refineFactor=16, int polyOrder=3,
- int derivative=0, bool storeCoefs=true);
- ~FftPolyInterpolator1D();
- virtual std::complex<double> F(double u);
- virtual std::complex<double> Dfdu(double u);
- virtual std::complex<double> Dfdfdudu(double u);
- std::complex<double> Integrate();
-};
-
-// FFT + spline interpolation on refined grid
-class FftSplineInterpolator1D : public FftPolyInterpolator1D {
- protected:
- // 1D spline interpolation
- std::complex<double> _SplineInterp(double *pnts,
- std::complex<double> *vals,
- std::complex<double> *deriv,
- double t);
- public:
- FftSplineInterpolator1D(const std::vector<double> &u,
- const std::vector< std::complex<double> > &data,
- int refineFactor=16);
- ~FftSplineInterpolator1D(){}
- virtual std::complex<double> F(double u);
-};
-
-#endif
diff --git a/contrib/FourierModel/Interpolator2D.cpp b/contrib/FourierModel/Interpolator2D.cpp
deleted file mode 100755
index e3fe9c6dad2864932e71bf6f9fa1f11cfc02f9c4..0000000000000000000000000000000000000000
--- a/contrib/FourierModel/Interpolator2D.cpp
+++ /dev/null
@@ -1,754 +0,0 @@
-#include <vector>
-#include <math.h>
-#include <complex>
-#include "Interpolator2D.h"
-
-extern "C" {
- void zgelss_(int &,int &,int &,std::complex<double> *,int &,
- std::complex<double> *,int &,double *,double &,int &,
- std::complex<double> *,int &,double *,int &);
-}
-
-FftPolyInterpolator2D::
-FftPolyInterpolator2D(const std::vector<double> &u,
- const std::vector<double> &v,
- const std::vector< std::vector< std::complex<double> > >
- &data, int refineFactor, int polyOrder, int derivative)
- : _refineFactor(refineFactor), _polyOrder(polyOrder), _derivative(derivative),
- _uFine(0), _vFine(0), _fineData(0), _fineDerivU(0), _fineDerivV(0),
- _fineDerivUU(0), _fineDerivVV(0), _fineDerivUV(0)
-{
- _uIntervals = u.size() - 1;
- _uFineIntervals = _refineFactor * _uIntervals;
- _vIntervals = v.size() - 1;
- _vFineIntervals = _refineFactor * _vIntervals;
-
- // min/max of input mesh
- if(u[_uIntervals] > u[0]) {
- _uMin = u[0];
- _uMax = u[_uIntervals];
- }
- else {
- _uMin = u[_uIntervals];
- _uMax = u[0];
- }
- if(v[_vIntervals] > u[0]) {
- _vMin = v[0];
- _vMax = v[_vIntervals];
- }
- else {
- _vMin = u[_vIntervals];
- _vMax = u[0];
- }
-
- // signed length of interval
- _uLength = u[_uIntervals] - u[0];
- _vLength = v[_vIntervals] - v[0];
-
- // Sanity checks
- if(std::abs(data[0][0] - data[_uIntervals][0]) > 1.e-10 ||
- std::abs(data[0][0] - data[0][_vIntervals]) > 1.e-10) {
- Msg::Warning("Data is not periodic");
- }
- if(fabs(u[1] - u[0]) - fabs(u[2] - u[1]) > 1.e-10 ||
- fabs(v[1] - v[0]) - fabs(v[2] - v[1]) > 1.e-10){
- Msg::Fatal("Input grid is not equispaced");
- }
- if((_uFineIntervals < _polyOrder) || (_vFineIntervals < _polyOrder)){
- Msg::Fatal("Number of intervals is smaller than polynomial order (%d < %d or %d < %d)",
- _uFineIntervals, _polyOrder, _vFineIntervals, _polyOrder);
- }
-
- // Compute the fine grid spacing (assuming the u and v are equally spaced)
- _hUFine = (u[1] - u[0]) / _refineFactor;
- _hVFine = (v[1] - v[0]) / _refineFactor;
-
- // Compute the points on the fine grid
- _uFine = new double[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++)
- _uFine[i] = u[0] + i * _hUFine;
-
- _vFine = new double[_vFineIntervals + 1];
- for(int i = 0; i < _vFineIntervals + 1; i++)
- _vFine[i] = v[0] + i * _hVFine;
-
- // Initialize _fineData to zero
- _fineData = new std::complex<double>*[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++){
- _fineData[i] = new std::complex<double>[_vFineIntervals + 1];
- for(int j = 0; j < _vFineIntervals + 1; j++)
- _fineData[i][j] = 0.;
- }
-
- std::complex<double> **forwardData = new std::complex<double>*[_uIntervals + 1];
- for(int i = 0; i < _uIntervals + 1; i++)
- forwardData[i] = new std::complex<double>[_vIntervals + 1];
-
- if (_refineFactor == 1){ // No Fourier interpolation
- // Copy data into _fineData
- for(int i = 0; i < _uFineIntervals + 1; i++)
- for(int j = 0; j < _vFineIntervals + 1; j++)
- _fineData[i][j] = data[i][j];
- }
- else {
- // Compute the forward FFT of the two dimensional data (only uses
- // data in i=0,...,uIntervals-1 and j=0,...,vIntervals-1 and
- // scales by 1/(uIntervals*vIntervals))
- for(int i = 0; i < _uIntervals + 1; i++)
- for(int j = 0; j < _vIntervals + 1; j++)
- forwardData[i][j] = data[i][j];
- _ForwardFft(_uIntervals + 1, _vIntervals + 1, forwardData);
-
- // Copy the Fourier coefficients into _fineData
- int halfIndexU = _uIntervals / 2;
- int halfIndexV = _vIntervals / 2;
- for(int i = 0; i < _uIntervals; i++)
- for(int j = 0; j < _vIntervals; j++){
- int iFine = i;
- int jFine = j;
- if(i > halfIndexU)
- iFine = _uFineIntervals + (i - _uIntervals);
- if(j > halfIndexV)
- jFine = _vFineIntervals + (j - _vIntervals);
- _fineData[iFine][jFine] = forwardData[i][j];
- }
-
- // Compute inverse FFT on _fineData to interpolate the data on the
- // fine grid
- _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineData);
- }
-
- // Check if we need to interpolate the derivative(s)
- if(_derivative){
- if(_refineFactor == 1){
- // No Fourier interpolation was performed, so we need to do it here...
- for(int i = 0; i < _uIntervals + 1; i++)
- for(int j = 0; j < _vIntervals + 1; j++)
- forwardData[i][j] = data[i][j];
- _ForwardFft(_uIntervals + 1, _vIntervals + 1, forwardData);
- }
-
- // Initialize _fineDeriv and _fineDeriv2 to zero
- if(_derivative & 1){
- _fineDerivU = new std::complex<double>*[_uFineIntervals + 1];
- _fineDerivV = new std::complex<double>*[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++){
- _fineDerivU[i] = new std::complex<double>[_vFineIntervals + 1];
- _fineDerivV[i] = new std::complex<double>[_vFineIntervals + 1];
- for(int j = 0; j < _vFineIntervals + 1; j++){
- _fineDerivU[i][j] = 0.;
- _fineDerivV[i][j] = 0.;
- }
- }
- }
-
- if(_derivative & 2){
- _fineDerivUU = new std::complex<double>*[_uFineIntervals + 1];
- _fineDerivVV = new std::complex<double>*[_uFineIntervals + 1];
- _fineDerivUV = new std::complex<double>*[_uFineIntervals + 1];
- for(int i = 0; i < _uFineIntervals + 1; i++){
- _fineDerivUU[i] = new std::complex<double>[_vFineIntervals + 1];
- _fineDerivVV[i] = new std::complex<double>[_vFineIntervals + 1];
- _fineDerivUV[i] = new std::complex<double>[_vFineIntervals + 1];
- for(int j = 0; j < _vFineIntervals + 1; j++){
- _fineDerivUU[i][j] = 0.;
- _fineDerivVV[i][j] = 0.;
- _fineDerivUV[i][j] = 0.;
- }
- }
- }
-
- // Copy the Fourier coefficients into _fineDeriv and _fineDeriv2
- std::complex<double> I(0., 1.);
- int halfIndexU = _uIntervals / 2;
- int halfIndexV = _vIntervals / 2;
- for(int i = 0; i < _uIntervals; i++){
- for(int j = 0; j < _vIntervals; j++){
- int iFine = i;
- int jFine = j;
- double kU = i;
- double kV = j;
- if(i > halfIndexU){
- kU = i - _uIntervals;
- iFine = _uFineIntervals + (i - _uIntervals);
- }
- if(j > halfIndexV){
- kV = j - _vIntervals;
- jFine = _vFineIntervals + (j - _vIntervals);
- }
- if(_derivative & 1){
- _fineDerivU[iFine][jFine] = (2. * M_PI * kU * I) * forwardData[i][j] / _uLength;
- _fineDerivV[iFine][jFine] = (2. * M_PI * kV * I) * forwardData[i][j] / _vLength;
- }
- if(_derivative & 2){
- _fineDerivUU[iFine][jFine] =
- -(4. * M_PI * M_PI * kU * kU) * forwardData[i][j] / (_uLength * _uLength);
- _fineDerivVV[iFine][jFine] =
- -(4. * M_PI * M_PI * kV * kV) * forwardData[i][j] / (_vLength * _vLength);
- _fineDerivUV[iFine][jFine] =
- -(4. * M_PI * M_PI * kU * kV) * forwardData[i][j] / (_uLength * _vLength);
- }
- }
- }
-
- // Perform an inverse FFT on _fineDeriv to interpolate the
- // derivative to the fine grid
- if(_derivative & 1){
- _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivU);
- _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivV);
- }
- if(_derivative & 2){
- _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivUU);
- _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivVV);
- _BackwardFft(_uFineIntervals + 1, _vFineIntervals + 1, _fineDerivUV);
- }
- }
-
- for(int i = 0; i < _uIntervals + 1; i++)
- delete [] forwardData[i];
- delete [] forwardData;
-
- // Initialize interpolation variables
- _tmpSpace = new double[_polyOrder + 1];
- _tmpValsReal = new double[_polyOrder + 1];
- _tmpValsImag = new double[_polyOrder + 1];
- _tmpInterpReal = new double[_polyOrder + 1];
- _tmpInterpImag = new double[_polyOrder + 1];
-}
-
-FftPolyInterpolator2D::~FftPolyInterpolator2D()
-{
- delete[] _uFine;
- delete[] _vFine;
-
- for(int i = 0; i < _uFineIntervals + 1; i++)
- delete [] _fineData[i];
- delete [] _fineData;
-
- if(_fineDerivU){
- for(int i = 0; i < _uFineIntervals + 1; i++)
- delete [] _fineDerivU[i];
- delete [] _fineDerivU;
- }
- if(_fineDerivV){
- for(int i = 0; i < _uFineIntervals + 1; i++)
- delete [] _fineDerivV[i];
- delete [] _fineDerivV;
- }
- if(_fineDerivUU){
- for(int i = 0; i < _uFineIntervals + 1; i++)
- delete [] _fineDerivUU[i];
- delete [] _fineDerivUU;
- }
- if(_fineDerivVV){
- for(int i = 0; i < _uFineIntervals + 1; i++)
- delete [] _fineDerivVV[i];
- delete [] _fineDerivVV;
- }
- if(_fineDerivUV){
- for(int i = 0; i < _uFineIntervals + 1; i++)
- delete [] _fineDerivUV[i];
- delete [] _fineDerivUV;
- }
-
- delete[] _tmpSpace;
- delete[] _tmpValsReal;
- delete[] _tmpValsImag;
- delete[] _tmpInterpReal;
- delete[] _tmpInterpImag;
-}
-
-int FftPolyInterpolator2D::_forwardSizeU = 0;
-int FftPolyInterpolator2D::_forwardSizeV = 0;
-int FftPolyInterpolator2D::_backwardSizeU = 0;
-int FftPolyInterpolator2D::_backwardSizeV = 0;
-fftw_plan FftPolyInterpolator2D::_forwardPlan;
-fftw_plan FftPolyInterpolator2D::_backwardPlan;
-fftw_complex *FftPolyInterpolator2D::_forwardData = 0;
-fftw_complex *FftPolyInterpolator2D::_backwardData = 0;
-
-void FftPolyInterpolator2D::_SetForwardPlan(int n, int m)
-{
- if(n != _forwardSizeU || m != _forwardSizeV){
- if(_forwardSizeU || _forwardSizeV){
- fftw_destroy_plan(_forwardPlan);
- fftw_free(_forwardData);
- }
- _forwardSizeU = n;
- _forwardSizeV = m;
- _forwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) *
- _forwardSizeU * _forwardSizeV);
- _forwardPlan = fftw_plan_dft_2d(_forwardSizeU, _forwardSizeV,
- _forwardData, _forwardData,
- FFTW_FORWARD, FFTW_ESTIMATE);
- }
-}
-
-void FftPolyInterpolator2D::_SetBackwardPlan(int n, int m)
-{
- if(n != _backwardSizeU || m != _backwardSizeV){
- if(_backwardSizeU || _backwardSizeV){
- fftw_destroy_plan(_backwardPlan);
- fftw_free(_backwardData);
- }
- _backwardSizeU = n;
- _backwardSizeV = m;
- _backwardData = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) *
- _backwardSizeU * _backwardSizeV);
- _backwardPlan = fftw_plan_dft_2d(_backwardSizeU, _backwardSizeV,
- _backwardData, _backwardData,
- FFTW_BACKWARD, FFTW_ESTIMATE);
- }
-}
-
-void FftPolyInterpolator2D::_ForwardFft(int n, int m, std::complex<double> **fftData)
-{
- // Initialize fftw plan and array (ignoring the last row and column
- // of fftData, which should just be the periodic extension)
- _SetForwardPlan(n - 1, m - 1);
- int k = 0;
- for(int i = 0; i < n - 1; i++){
- for(int j = 0; j < m - 1; j++){
- _forwardData[k][0] = fftData[i][j].real();
- _forwardData[k][1] = fftData[i][j].imag();
- k++;
- }
- }
-
- // Perform forward FFT
- fftw_execute(_forwardPlan);
-
- // Copy data back into fftData and scale by 1/((n-1) * (m-1))
- double s = 1. / (double)((n - 1) * (m - 1));
- k = 0;
- for(int i = 0; i < n - 1; i++){
- for(int j = 0; j < m - 1; j++){
- fftData[i][j] = s * std::complex<double>(_forwardData[k][0], _forwardData[k][1]);
- k++;
- }
- }
-}
-
-void FftPolyInterpolator2D::_BackwardFft(int n, int m, std::complex<double> **fftData)
-{
- // Initialize fftw plan and array (ignoring last row and column of fftData)
- _SetBackwardPlan(n - 1, m - 1);
- int k = 0;
- for(int i = 0; i < n - 1; i++){
- for(int j = 0; j < m - 1; j++){
- _backwardData[k][0] = fftData[i][j].real();
- _backwardData[k][1] = fftData[i][j].imag();
- k++;
- }
- }
-
- // Perform backward FFT
- fftw_execute(_backwardPlan);
-
- // Copy data back into fftData
- k = 0;
- for(int i = 0; i < n - 1; i++){
- for(int j = 0; j < m - 1; j++){
- fftData[i][j] = std::complex<double>(_backwardData[k][0], _backwardData[k][1]);
- k++;
- }
- }
-
- // Fill in last row and column with copy of first row and column
- for(int i = 0; i <= n - 1; i++)
- fftData[i][m - 1] = fftData[i][0];
- for(int j = 0; j <= m - 1; j++)
- fftData[n - 1][j] = fftData[0][j];
-}
-
-double FftPolyInterpolator2D::_PolyInterp(double start, double h, const double *vals,
- double t, double *space)
-{
- // Neville's Algorithm from Stoer and Bulirsch, Section 2.1.2
- double tmp = (t - start)/h;
- for(int i = 0; i <= _polyOrder; i++){
- space[i] = vals[i];
- for(int j = i - 1; j >= 0; j--)
- space[j] = space[j + 1] + (space[j + 1] - space[j]) * ((tmp - i) / (i - j));
- }
- return space[0];
-}
-
-std::complex<double> FftPolyInterpolator2D::_Interpolate(double u, double v,
- int uDer, int vDer)
-{
- if( ( (uDer == 2 || vDer == 2 || (uDer == 1 && vDer == 1)) && !(_derivative & 2) ) ||
- ( (uDer == 1 || vDer == 1) && !(_derivative & 1) ) ||
- (uDer < 0 || uDer > 2 || vDer < 0 || vDer > 2) ){
- Msg::Error("Derivative data not available: check contructor call %d %d %d",uDer,vDer,_derivative);
- return 0.;
- }
-
- double epsilon = 1e-12;
- if(u < _uMin - epsilon || u > _uMax + epsilon){
- Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
- "not in [%g,%g]x[%g,%g]", u, v, _uMin, _uMax, _vMin, _vMax);
- return 0.;
- }
-
- // Choose uIndexStart and vIndexStart so that (u,v) is centered in
- // the polynomial interpolation grid
- int uIndexStart = (int)floor((u - _uFine[0]) / _hUFine) - _polyOrder / 2;
- if(uIndexStart < 0)
- uIndexStart = 0;
- else if((uIndexStart + _polyOrder) > _uFineIntervals)
- uIndexStart = _uFineIntervals - _polyOrder;
-
- int vIndexStart = (int)floor((v - _vFine[0]) / _hVFine) - _polyOrder / 2;
- if(vIndexStart < 0)
- vIndexStart = 0;
- else if((vIndexStart + _polyOrder) > _vFineIntervals)
- vIndexStart = _vFineIntervals - _polyOrder;
-
- // Interpolate to find value at (u,v)
- double start = _vFine[vIndexStart];
- for(int i = 0; i <= _polyOrder; i++){
- for(int j = 0; j <= _polyOrder; j++){
- std::complex<double> tmp;
- if(uDer == 0 && vDer == 0)
- tmp = _fineData[uIndexStart + i][vIndexStart + j];
- else if(uDer == 1 && vDer == 0)
- tmp = _fineDerivU[uIndexStart + i][vIndexStart + j];
- else if(uDer == 0 && vDer == 1)
- tmp = _fineDerivV[uIndexStart + i][vIndexStart + j];
- else if(uDer == 2 && vDer == 0)
- tmp = _fineDerivUU[uIndexStart + i][vIndexStart + j];
- else if(uDer == 0 && vDer == 2)
- tmp = _fineDerivVV[uIndexStart + i][vIndexStart + j];
- else
- tmp = _fineDerivUV[uIndexStart + i][vIndexStart + j];
- _tmpValsReal[j] = tmp.real();
- _tmpValsImag[j] = tmp.imag();
- }
- _tmpInterpReal[i] = _PolyInterp(start, _hVFine, _tmpValsReal, v, _tmpSpace);
- _tmpInterpImag[i] = _PolyInterp(start, _hVFine, _tmpValsImag, v, _tmpSpace);
- }
-
- start = _uFine[uIndexStart];
- return std::complex<double>(_PolyInterp(start, _hUFine, _tmpInterpReal, u, _tmpSpace),
- _PolyInterp(start, _hUFine, _tmpInterpImag, u, _tmpSpace));
-}
-
-std::complex<double> FftPolyInterpolator2D::F(double u, double v)
-{
- return _Interpolate(u, v, 0, 0);
-}
-
-std::complex<double> FftPolyInterpolator2D::Dfdu(double u, double v)
-{
- return _Interpolate(u, v, 1, 0);
-}
-
-std::complex<double> FftPolyInterpolator2D::Dfdv(double u, double v)
-{
- return _Interpolate(u, v, 0, 1);
-}
-
-std::complex<double> FftPolyInterpolator2D::Dfdfdudu(double u, double v)
-{
- return _Interpolate(u, v, 2, 0);
-}
-
-std::complex<double> FftPolyInterpolator2D::Dfdfdvdv(double u, double v)
-{
- return _Interpolate(u, v, 0, 2);
-}
-
-std::complex<double> FftPolyInterpolator2D::Dfdfdudv(double u, double v)
-{
- return _Interpolate(u, v, 1, 1);
-}
-
-
-
-FourierContinuationInterpolator2D::FourierContinuationInterpolator2D
-(const std::vector<double> &u, const std::vector<double> &v,
- const std::vector< std::complex<double> > &data,
- int derivative, int uModes, int vModes, double periodU, double periodV)
- : _derivative(derivative), _uModes(uModes), _vModes(vModes),
- _periodU(periodU), _periodV(periodV),
- _coeffData(0), _coeffDerivU(0), _coeffDerivV(0), _coeffDerivUU(0),
- _coeffDerivVV(0), _coeffDerivUV(0)
-{
- // sanity check
- if (!((u.size() == v.size()) && (v.size() == data.size()))) {
- Msg::Fatal("Input data of unequal size.");
- }
- _nData = data.size();
- _nModes = _uModes * _vModes;
-
- for (int i=0;i<_nData;i++) {
- _u.push_back(u[i]);
- _v.push_back(v[i]);
- }
-
- // Initialize _Data to zero
- _coeffData = new std::complex<double>*[_uModes];
- for(int j = 0; j < _uModes; j++){
- _coeffData[j] = new std::complex<double>[_vModes];
- for(int k = 0; k < _vModes; k++)
- _coeffData[j][k] = 0.;
- }
-
- if ((_uModes % 2) == 0) {
- _uModesLower = - _uModes/2;
- _uModesUpper = _uModes/2 - 1;
- }
- else {
- _uModesLower = - (_uModes - 1)/2;
- _uModesUpper = (_uModes - 1)/2;
- }
- if ((_vModes % 2) == 0) {
- _vModesLower = - _vModes/2;
- _vModesUpper = _vModes/2 - 1;
- }
- else {
- _vModesLower = - (_vModes - 1)/2;
- _vModesUpper = (_vModes - 1)/2;
- }
-
- std::complex<double> *LSRhs = new std::complex<double> [_nData];
- std::complex<double> *LSMatrix = new std::complex<double> [_nModes * _nData];
-
- for (int i=0;i<_nData;i++)
- LSRhs[i] = data[i];
-
- for (int j=0;j<_uModes;j++)
- for (int k=0;k<_vModes;k++)
- for (int i=0;i<_nData;i++)
- LSMatrix[i+_nData*(k+_vModes*j)] = std::complex<double>
- (cos(2*M_PI * ((j + _uModesLower) * u[i] / _periodU +
- (k + _vModesLower) * v[i] / _periodV)),
- sin(2*M_PI * ((j + _uModesLower) * u[i] / _periodU +
- (k + _vModesLower) * v[i] / _periodV)));
-
- // some parameters for the lease square solvers "zgelss"
- int info, rank;
- int lwork = 66*_nData, one = 1;
- double rcond = -1.;
- double *s = new double [_nModes];
- double *rwork = new double [5*_nModes];
- std::complex<double> *work = new std::complex<double> [lwork];
-
- zgelss_(_nData,_nModes,one,&LSMatrix[0],_nData,&LSRhs[0],_nData,&s[0],rcond,
- rank,&work[0],lwork,&rwork[0],info);
-
- for(int j = 0; j < _uModes; j++)
- for(int k = 0; k < _vModes; k++)
- _coeffData[j][k] = LSRhs[k+_vModes*j];
-
- // deleting lease square arrays
- delete[] s, rwork, work;
- delete[] LSMatrix, LSRhs;
-
- // Check if we need to interpolate the derivative(s)
- if(_derivative){
- // Initialize _fineDeriv and _fineDeriv2 to zero
- if(_derivative & 1){
- _coeffDerivU = new std::complex<double>*[_uModes];
- _coeffDerivV = new std::complex<double>*[_uModes];
- for(int j = 0; j < _uModes; j++){
- _coeffDerivU[j] = new std::complex<double>[_vModes];
- _coeffDerivV[j] = new std::complex<double>[_vModes];
- for(int k = 0; k < _vModes; k++){
- _coeffDerivU[j][k] = 0.;
- _coeffDerivV[j][k] = 0.;
- }
- }
- }
-
- if(_derivative & 2){
- _coeffDerivUU = new std::complex<double>*[_uModes];
- _coeffDerivVV = new std::complex<double>*[_uModes];
- _coeffDerivUV = new std::complex<double>*[_uModes];
- for(int j = 0; j < _uModes; j++){
- _coeffDerivUU[j] = new std::complex<double>[_vModes];
- _coeffDerivVV[j] = new std::complex<double>[_vModes];
- _coeffDerivUV[j] = new std::complex<double>[_vModes];
- for(int k = 0; k < _vModes; k++){
- _coeffDerivUU[j][k] = 0.;
- _coeffDerivVV[j][k] = 0.;
- _coeffDerivUV[j][k] = 0.;
- }
- }
- }
-
- // Copy the Fourier coefficients into _coeffDeriv and _coeffDeriv2
- std::complex<double> I(0., 1.);
- for(int j = 0; j < _uModes; j++){
- for(int k = 0; k < _vModes; k++){
- int J = j+_uModesLower;
- int K = k+_vModesLower;
- if(_derivative & 1){
- _coeffDerivU[j][k] = (2 * M_PI * J * I / _periodU) *_coeffData[j][k];
- _coeffDerivV[j][k] = (2 * M_PI * K * I / _periodV) *_coeffData[j][k];
- }
- if(_derivative & 2){
- _coeffDerivUU[j][k] = - (4 * M_PI * M_PI * J * J /
- (_periodU * _periodU)) * _coeffData[j][k];
- _coeffDerivVV[j][k] = - (4 * M_PI * M_PI * K * K /
- (_periodV * _periodV)) * _coeffData[j][k];
- _coeffDerivUV[j][k] = - (4 * M_PI * M_PI * J * K /
- (_periodU * _periodV)) * _coeffData[j][k];
- }
- }
- }
- }
- // Initialize interpolation variables
- _tmpCoeff = std::vector< std::complex<double> >(_vModes);
- _tmpInterp = std::vector< std::complex<double> >(_uModes);
-}
-
-FourierContinuationInterpolator2D::~FourierContinuationInterpolator2D()
-{
- for(int j = 0; j < _uModes; j++)
- delete [] _coeffData[j];
- delete [] _coeffData;
-
- if(_coeffDerivU){
- for(int j = 0; j < _uModes; j++)
- delete [] _coeffDerivU[j];
- delete [] _coeffDerivU;
- }
- if(_coeffDerivV){
- for(int j = 0; j < _uModes; j++)
- delete [] _coeffDerivV[j];
- delete [] _coeffDerivV;
- }
- if(_coeffDerivUU){
- for(int j = 0; j < _uModes; j++)
- delete [] _coeffDerivUU[j];
- delete [] _coeffDerivUU;
- }
- if(_coeffDerivVV){
- for(int j = 0; j < _uModes; j++)
- delete [] _coeffDerivVV[j];
- delete [] _coeffDerivVV;
- }
- if(_coeffDerivUV){
- for(int j = 0; j < _uModes; j++)
- delete [] _coeffDerivUV[j];
- delete [] _coeffDerivUV;
- }
-}
-
-std::complex<double> FourierContinuationInterpolator2D::
-_PolyEval(std::vector< std::complex<double> > _coeff, std::complex<double> x)
-{
- int _polyOrder = _coeff.size()-1;
- std::complex<double> out = 0.;
-
- out = x * _coeff[_polyOrder];
- for (int i = _polyOrder - 1; i > 0; i--)
- out = x * (out + _coeff[i]);
- out = out + _coeff[0];
-
- return out;
-}
-
-std::complex<double> FourierContinuationInterpolator2D::
- _Interpolate(double u, double v, int uDer, int vDer)
-{
- //Msg::Info("%d %d %d",uDer,vDer,_derivative);
- if (((uDer==2 || vDer==2 || (uDer==1 && vDer==1)) && !(_derivative & 2) ) ||
- ((uDer==1 || vDer==1) && !(_derivative & 1)) ||
- (uDer<0 || uDer>2 || vDer<0 || vDer>2) ) {
- Msg::Error("Derivative data not available: check contructor call %d %d %d",
- uDer,vDer,_derivative);
- return 0.;
- }
-
- double epsilon = 1e-12;
- if(u < 0. - epsilon || u > 1. + epsilon){
- Msg::Error("Trying to interpolate outside interval: (u,v)=(%.16g,%.16g) "
- "not in [%g,%g]x[%g,%g]", u, v, 0., 1., 0., 1.);
- return 0.;
- }
-
- // Interpolate to find value at (u,v)
- for(int j = 0; j < _uModes; j++){
- for(int k = 0; k < _vModes; k++){
- std::complex<double> tmp;
- if(uDer == 0 && vDer == 0)
- tmp = _coeffData[j][k];
- else if(uDer == 1 && vDer == 0)
- tmp = _coeffDerivU[j][k];
- else if(uDer == 0 && vDer == 1)
- tmp = _coeffDerivV[j][k];
- else if(uDer == 2 && vDer == 0)
- tmp = _coeffDerivUU[j][k];
- else if(uDer == 0 && vDer == 2)
- tmp = _coeffDerivVV[j][k];
- else
- tmp = _coeffDerivUV[j][k];
- _tmpCoeff[k] = tmp;
- }
- std::complex<double> y(cos(2 * M_PI * v / _periodV),
- sin(2 * M_PI * v / _periodV));
- _tmpInterp[j] = _PolyEval(_tmpCoeff, y);
- _tmpInterp[j] *= std::complex<double>
- (cos(2 * M_PI * _vModesLower * v / _periodV),
- sin(2 * M_PI * _vModesLower * v / _periodV));
- }
- std::complex<double> x(cos(2 * M_PI * u / _periodU),
- sin(2 * M_PI * u / _periodU));
- return _PolyEval(_tmpInterp, x) * std::complex<double>
- (cos(2 * M_PI * _uModesLower * u / _periodU),
- sin(2 * M_PI * _uModesLower * u / _periodU));
-}
-
-std::complex<double> FourierContinuationInterpolator2D::F(double u, double v)
-{
- return _Interpolate(u, v, 0, 0);
-}
-
-std::complex<double> FourierContinuationInterpolator2D::Dfdu(double u, double v)
-{
- return _Interpolate(u, v, 1, 0);
-}
-
-std::complex<double> FourierContinuationInterpolator2D::Dfdv(double u, double v)
-{
- return _Interpolate(u, v, 0, 1);
-}
-
-std::complex<double> FourierContinuationInterpolator2D::Dfdfdudu(double u, double v)
-{
- return _Interpolate(u, v, 2, 0);
-}
-
-std::complex<double> FourierContinuationInterpolator2D::Dfdfdvdv(double u, double v)
-{
- return _Interpolate(u, v, 0, 2);
-}
-
-std::complex<double> FourierContinuationInterpolator2D::Dfdfdudv
-(double u, double v)
-{
- return _Interpolate(u, v, 1, 1);
-}
-
-bool FourierContinuationInterpolator2D::IsPointInInterpolationRange
-(double u, double v)
-{
- double nbdRadius = 1.e-2;
- bool status = false;
- for (int i=0;i<_nData;i++) {
- double tmp = sqrt((u - _u[i])*(u - _u[i])+(v - _v[i])*(v - _v[i]));
- if (tmp < nbdRadius) {
- status = true;
- break;
- }
- }
- return status;
-}
-
-int FourierContinuationInterpolator2D::GetDataSize()
-{
- return _nData;
-}
diff --git a/contrib/FourierModel/Interpolator2D.h b/contrib/FourierModel/Interpolator2D.h
deleted file mode 100755
index 1b55eb5fc156b74dabbda72a383af2967f1b88cd..0000000000000000000000000000000000000000
--- a/contrib/FourierModel/Interpolator2D.h
+++ /dev/null
@@ -1,166 +0,0 @@
-#ifndef _INTERPOLATOR_2D_H_
-#define _INTERPOLATOR_2D_H_
-
-#include <vector>
-#include <complex>
-#include <fftw3.h>
-#include "Message.h"
-
-// The base class for the 2D interpolators
-class Interpolator2D {
- public:
- Interpolator2D(){}
- virtual ~Interpolator2D(){}
- // interpolates the data at u,v
- virtual std::complex<double> F(double u, double v) = 0;
- // interpolates the first u derivative of the data at u,v
- virtual std::complex<double> Dfdu(double u, double v)
- {
- Msg::Error("First derivative not implemented for this interpolator");
- return 0.;
- }
- // interpolates the first v derivative of the data at u,v
- virtual std::complex<double> Dfdv(double u, double v)
- {
- Msg::Error("First derivative not implemented for this interpolator");
- return 0.;
- }
- // interpolates the second u derivative of the data at u,v
- virtual std::complex<double> Dfdfdudu(double u, double v)
- {
- Msg::Error("Second derivative not implemented for this interpolator");
- return 0.;
- }
- // interpolates the second v derivative of the data at u,v
- virtual std::complex<double> Dfdfdvdv(double u, double v)
- {
- Msg::Error("Second derivative not implemented for this interpolator");
- return 0.;
- }
- // interpolates the second cross derivative of the data at u,v
- virtual std::complex<double> Dfdfdudv(double u, double v)
- {
- Msg::Error("Second derivative not implemented for this interpolator");
- return 0.;
- }
- // checks if the interpolation point is good enough
- virtual bool IsPointInInterpolationRange(double u, double v)
- {
- Msg::Error("Goodness check not implemented for this interpolator");
- return 0;
- }
- // returns the size of interpolation data
- virtual int GetDataSize()
- {
- Msg::Error("GetDataSize not implemented for this interpolator");
- return 0;
- }
-};
-
-// FFT + polynomial interpolation on refined grid (assumes that the
-// input grid is equally spaced and the input data is periodic)
-class FftPolyInterpolator2D : public Interpolator2D {
- private:
- // refinement factor (e.g. 16; if equal to 1, we just perform
- // standard piecewise polynomial interpolation, without any FFTs)
- int _refineFactor;
- // order of the interpolating polynomial
- int _polyOrder;
- // signed length of the interval
- double _uLength, _vLength;
- // temporary interpolation variables
- double *_tmpSpace, *_tmpValsReal, *_tmpValsImag, *_tmpInterpReal, *_tmpInterpImag;
- // 1D polynomial interpolation routine
- double _PolyInterp(double start, double h, const double *vals,
- double t, double *space);
- // interpolation wrapper
- std::complex<double> _Interpolate(double u, double v, int uDer=0, int vDer=0);
- // persistent fftw data (used to avoid recomputing the FFTW plans
- // and reallocating memory every time)
- static int _forwardSizeU, _forwardSizeV;
- static int _backwardSizeU, _backwardSizeV;
- static fftw_plan _forwardPlan, _backwardPlan;
- static fftw_complex *_forwardData, *_backwardData;
- void _SetForwardPlan(int n, int m);
- void _SetBackwardPlan(int n, int m);
- protected:
- // number of intervals in the original and in the refined mesh
- int _uIntervals, _vIntervals, _uFineIntervals, _vFineIntervals;
- // fine mesh spacing
- double _hUFine, _hVFine;
- // min/max of input mesh
- double _uMin, _uMax, _vMin, _vMax;
- // bitfield telling if we also interpolate the derivative(s)
- int _derivative;
- // grid points of the refined mesh
- double *_uFine, *_vFine;
- // data (and its first 2 derivatives) on the refined regular grid
- std::complex<double> **_fineData, **_fineDerivU, **_fineDerivV;
- std::complex<double> **_fineDerivUU, **_fineDerivVV, **_fineDerivUV;
- // This routine computes the forward FFT (ignoring the last row and
- // column of the data)
- void _ForwardFft(int n, int m, std::complex<double> **fftData);
- // This routine computes the inverse FFT (ignoring the last row and
- // column in the array), returns values in entire array (by using
- // the periodic extension)
- void _BackwardFft(int n, int m, std::complex<double> **fftData);
- public:
- FftPolyInterpolator2D(const std::vector<double> &u,
- const std::vector<double> &v,
- const std::vector<std::vector<std::complex<double> > >
- &data, int refineFactor=16, int polyOrder=3,
- int derivative=0);
- ~FftPolyInterpolator2D();
- virtual std::complex<double> F(double u, double v);
- virtual std::complex<double> Dfdu(double u, double v);
- virtual std::complex<double> Dfdv(double u, double v);
- virtual std::complex<double> Dfdfdudu(double u, double v);
- virtual std::complex<double> Dfdfdvdv(double u, double v);
- virtual std::complex<double> Dfdfdudv(double u, double v);
-};
-
-// Fourier Continuation interpolation
-class FourierContinuationInterpolator2D : public Interpolator2D {
- private:
- friend class Body;
- // u and v data
- std::vector<double> _u, _v;
- // temporary interpolation variables
- std::vector< std::complex<double> > _tmpCoeff, _tmpInterp;
- protected:
- // bitfield telling if we also interpolate the derivative(s)
- int _derivative;
- // Data Size
- int _nData;
- // Number of Fourier Modes
- int _uModes, _vModes, _nModes;
- // Period Information
- double _periodU, _periodV;
- // Limits of Fourier Series
- int _uModesLower, _uModesUpper, _vModesLower, _vModesUpper;
- // data (and its first 2 derivatives)
- std::complex<double> **_coeffData, **_coeffDerivU, **_coeffDerivV;
- std::complex<double> **_coeffDerivUU, **_coeffDerivVV, **_coeffDerivUV;
- // polynomial evaluator
- std::complex<double> _PolyEval(std::vector< std::complex<double> > _coeff,
- std::complex<double> x);
- // interpolation wrapper
- std::complex<double> _Interpolate(double u,double v,int uDer=0,int vDer=0);
- public:
- FourierContinuationInterpolator2D
- (const std::vector<double> &u, const std::vector<double> &v,
- const std::vector< std::complex<double> > &data,
- int derivative = 0, int uModes = 1, int vModes = 1,
- double periodU = 2, double periodV = 2);
- ~FourierContinuationInterpolator2D();
- virtual std::complex<double> F(double u, double v);
- virtual std::complex<double> Dfdu(double u, double v);
- virtual std::complex<double> Dfdv(double u, double v);
- virtual std::complex<double> Dfdfdudu(double u, double v);
- virtual std::complex<double> Dfdfdvdv(double u, double v);
- virtual std::complex<double> Dfdfdudv(double u, double v);
- virtual bool IsPointInInterpolationRange(double u, double v);
- virtual int GetDataSize();
-};
-
-#endif
diff --git a/contrib/FourierModel/Makefile b/contrib/FourierModel/Makefile
index c3f73e570b0f33ce4bd2763e5b308e690203e901..64704307c0983e7131f72e58e2238ed9fa50b873 100644
--- a/contrib/FourierModel/Makefile
+++ b/contrib/FourierModel/Makefile
@@ -75,9 +75,9 @@ FM_Face.o: FM_Face.cpp FM_Face.h Patch.h ProjectionSurface.h FM_Edge.h \
FM_Info.o: FM_Info.cpp FM_Info.h
FM_Reader.o: FM_Reader.cpp Message.h FM_Reader.h Curve.h \
IntersectionCurve.h Patch.h ProjectionSurface.h FM_Info.h ExactPatch.h \
- ContinuationPatch.h PartitionOfUnity.h \
- CylindricalProjectionSurface.h RevolvedParabolaProjectionSurface.h \
- Utils.h FM_Face.h FM_Edge.h FM_Vertex.h BlendOperator.h BlendedPatch.h
+ ContinuationPatch.h PartitionOfUnity.h CylindricalProjectionSurface.h \
+ RevolvedParabolaProjectionSurface.h Utils.h FM_Face.h FM_Edge.h \
+ FM_Vertex.h BlendOperator.h BlendedPatch.h
FM_Vertex.o: FM_Vertex.cpp
Message.o: Message.cpp Message.h
Utils.o: Utils.cpp Utils.h Message.h