36 #define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr 37 #define TWOPI 6.28318530717959 51 unsigned int dimCount = nn[idim];
52 if (!dimCount || (dimCount & (dimCount - 1)))
56 "fft::transform(complexField&, const labelList&, " 58 ) <<
"number of elements in direction " << idim
59 <<
" is not a power of 2" <<
endl 60 <<
" Number of elements in each direction = " << nn
67 label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
68 label ibit, k1, k2,
n, nprev, nrem, idim;
70 scalar theta, wi, wpi, wpr, wr, wtemp;
71 scalar*
data =
reinterpret_cast<scalar*
>(field.
begin()) - 1;
76 if (isign == REVERSE_TRANSFORM)
91 for (idim=ndim; idim>=1; idim--)
94 nrem = ntot/(n*nprev);
100 for (i2=1; i2<=ip2; i2+=ip1)
104 for (i1=i2; i1<=i2 + ip1 - 2; i1+=2)
106 for (i3=i1; i3<=ip3; i3+=ip2)
108 i3rev = i2rev + i3 - i2;
109 SWAP(data[i3], data[i3rev]);
110 SWAP(data[i3 + 1], data[i3rev + 1]);
116 while (ibit >= ip1 && i2rev > ibit)
130 theta = isign*
TWOPI/(ifp2/ip1);
131 wtemp =
sin(0.5*theta);
132 wpr = -2.0*wtemp*wtemp;
137 for (i3 = 1; i3 <= ifp1; i3 += ip1)
139 for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
141 for (i2 = i1; i2 <= ip3; i2 += ifp2)
145 tempr = scalar(wr*data[k2]) - scalar(wi*data[k2 + 1]);
146 tempi = scalar(wr*data[k2 + 1]) + scalar(wi*data[k2]);
147 data[k2] = data[k1] - tempr;
148 data[k2 + 1] = data[k1 + 1] - tempi;
150 data[k1 + 1] += tempi;
154 wr = (wtemp = wr)*wpr - wi*wpi + wr;
155 wi = wi*wpr + wtemp*wpi + wi;
165 if (isign == FORWARD_TRANSFORM)
173 scalar recRootN = 1.0/
sqrt(scalar(ntot));
177 field[i] *= recRootN;
197 transform(tfftField(), nn, FORWARD_TRANSFORM);
213 transform(tifftField(), nn, REVERSE_TRANSFORM);
237 tfftVectorField().replace
240 forwardTransform(tfield().
component(cmpt), nn)
246 return tfftVectorField;
266 tifftVectorField().replace
269 reverseTransform(tfield().
component(cmpt), nn)
275 return tifftVectorField;
dimensionedScalar sqrt(const dimensionedScalar &ds)
void clear() const
If object pointer points to valid object:
Multi-dimensional renumbering used in the Numerical Recipes fft routine.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void size(const label)
Override size to be inconsistent with allocated storage.
Database for solution data, solver performance and other reduced data.
Field< complex > complexField
Specialisation of Field<T> for complex.
static tmp< complexField > reverseTransform(const tmp< complexField > &field, const labelList &nn)
Number of components in this vector space.
iterator begin()
Return an iterator to begin traversing the UList.
Ostream & endl(Ostream &os)
Add newline and flush stream.
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
dimensionSet transform(const dimensionSet &)
void fftRenumber(List< complex > &data, const labelList &nn)
Pre-declare SubField and related Field type.
errorManip< error > abort(error &err)
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
A class for managing temporary objects.
static void transform(complexField &field, const labelList &nn, transformDirection fftDirection)
static tmp< complexField > forwardTransform(const tmp< complexField > &field, const labelList &nn)
dimensionedScalar sin(const dimensionedScalar &ds)