diff --git a/src/randomProcesses/fft/fft.C b/src/randomProcesses/fft/fft.C index fd44c0109df10e7ab98a7df1456789ee5d56b0d8..608fd2fd3ab52f894fafecc5a853d3d42828e42e 100644 --- a/src/randomProcesses/fft/fft.C +++ b/src/randomProcesses/fft/fft.C @@ -34,6 +34,56 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field) +{ + // Copy of input field for use by fftw + // - fftw requires non-const access to input and output... + scalarField fftInOut(field); + + const label n = fftInOut.size(); + const label nBy2 = n/2; + + // Using real to half-complex fftw 'kind' + fftw_plan plan = fftw_plan_r2r_1d + ( + n, + fftInOut.data(), + fftInOut.data(), + FFTW_R2HC, + FFTW_ESTIMATE + ); + + fftw_execute(plan); + + fftw_destroy_plan(plan); + + // field[0] = DC component + tmp<complexField> tresult(new complexField(nBy2 + 1)); + complexField& result = tresult.ref(); + + result[0].Re() = fftInOut[0]; + result[nBy2].Re() = fftInOut[nBy2]; + for (label i = 1; i < nBy2; ++i) + { + result[i].Re() = fftInOut[i]; + result[i].Im() = fftInOut[n - i]; + } + + return tresult; +} + + +Foam::tmp<Foam::complexField> fft::realTransform1D +( + const tmp<scalarField>& tfield +) +{ + tmp<complexField> tresult = realTransform1D(tfield()); + tfield.clear(); + return tresult; +} + + void fft::transform ( complexField& field, diff --git a/src/randomProcesses/fft/fft.H b/src/randomProcesses/fft/fft.H index 11cc07575e99b99c41f8942f31647223982d0360..7b62e44757e5414fa06c4bf76f82701b395e059e 100644 --- a/src/randomProcesses/fft/fft.H +++ b/src/randomProcesses/fft/fft.H @@ -44,7 +44,6 @@ SourceFiles #define fft_H #include "complexFields.H" -#include "UList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -63,6 +62,19 @@ public: }; + //- Transform real-value data + // - uses the fftw real to half-complex method + // - result size is field.size()/2 + 1 + static tmp<complexField> realTransform1D(const scalarField& field); + + + //- Transform real-value data + // - uses the fftw real to half-complex method + // - result size is field.size()/2 + 1 + static tmp<complexField> realTransform1D(const tmp<scalarField>& field); + + + //- Transform compex-value data static void transform ( complexField& field, diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C index e82d25b90fbef1fc4f72a70ea37adf37a1c4a6ed..57c46d7c0c7023cc7f5d41eb3e57ec83bb62d98d 100644 --- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C +++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C @@ -26,7 +26,6 @@ License #include "noiseFFT.H" #include "IFstream.H" #include "DynamicList.H" -#include "SubField.H" #include "mathematicalConstants.H" #include "HashSet.H" #include "fft.H" @@ -240,33 +239,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf const tmp<scalarField>& tpn ) const { - // Calculate the 2-sided fft - // Note: result not scaled - tmp<scalarField> tPn2 - ( - mag - ( - fft::reverseTransform - ( - ReComplexField(tpn), - List<int>(1, tpn().size()) - ) - ) - ); - - tpn.clear(); - - // Only storing the positive half - // Note: storing (N/2) values, DC component at position (0) - tmp<scalarField> tPn - ( - new scalarField - ( - scalarField::subField(tPn2(), tPn2().size()/2 + 1) - ) - ); - - return tPn; + return mag(fft::realTransform1D(tpn)); } diff --git a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C index 9d3e478c043cf52eeb10ec88316371475c5771cd..90f9971381fb767673c6ffacb9dd2bfbcf026a18 100644 --- a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C +++ b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C @@ -23,6 +23,10 @@ License \*---------------------------------------------------------------------------*/ +#include "SubField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + template<class Type> Foam::tmp<Foam::Field<Type>> Foam::windowModel::apply (