src/randomProcesses/fft/fft.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "fft.H"
  27 #include "fftRenumber.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 #define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr
  37 #define TWOPI 6.28318530717959
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 void fft::transform
  42 (
  43     complexField& field,
  44     const labelList& nn,
  45     transformDirection isign
  46 )
  47 {
  48     forAll(nn, idim)
  49     {
  50         // Check for power of two
  51         unsigned int dimCount = nn[idim];
  52         if (!dimCount || (dimCount & (dimCount - 1)))
  53         {
  54             FatalErrorIn
  55             (
  56                  "fft::transform(complexField&, const labelList&, "
  57                  "transformDirection)"
  58             )   << "number of elements in direction " << idim
  59                 << " is not a power of 2" << endl
  60                 << "    Number of elements in each direction = " << nn
  61                 << abort(FatalError);
  62         }
  63     }
  64
  65     const label ndim = nn.size();
  66
  67     label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
  68     label ibit, k1, k2, n, nprev, nrem, idim;
  69     scalar tempi, tempr;
  70     scalar theta, wi, wpi, wpr, wr, wtemp;
  71     scalar* data = reinterpret_cast<scalar*>(field.begin()) - 1;
  72
  73
  74     // if inverse transform : renumber before transform
  75
  76     if (isign == REVERSE_TRANSFORM)
  77     {
  78         fftRenumber(field, nn);
  79     }
  80
  81
  82     label ntot = 1;
  83     forAll(nn, idim)
  84     {
  85         ntot *= nn[idim];
  86     }
  87
  88
  89     nprev = 1;
  90
  91     for (idim=ndim; idim>=1; idim--)
  92     {
  93         n = nn[idim-1];
  94         nrem = ntot/(n*nprev);
  95         ip1 = nprev << 1;
  96         ip2 = ip1*n;
  97         ip3 = ip2*nrem;
  98         i2rev = 1;
  99
 100         for (i2=1; i2<=ip2; i2+=ip1)
 101         {
 102             if (i2 < i2rev)
 103             {
 104                 for (i1=i2; i1<=i2 + ip1 - 2; i1+=2)
 105                 {
 106                     for (i3=i1; i3<=ip3; i3+=ip2)
 107                     {
 108                         i3rev = i2rev + i3 - i2;
 109                         SWAP(data[i3], data[i3rev]);
 110                         SWAP(data[i3 + 1], data[i3rev + 1]);
 111                     }
 112                 }
 113             }
 114
 115             ibit = ip2 >> 1;
 116             while (ibit >= ip1 && i2rev > ibit)
 117             {
 118                 i2rev -= ibit;
 119                 ibit >>= 1;
 120             }
 121
 122             i2rev += ibit;
 123         }
 124
 125         ifp1 = ip1;
 126
 127         while (ifp1 < ip2)
 128         {
 129             ifp2 = ifp1 << 1;
 130             theta = isign*TWOPI/(ifp2/ip1);
 131             wtemp = sin(0.5*theta);
 132             wpr = -2.0*wtemp*wtemp;
 133             wpi = sin(theta);
 134             wr = 1.0;
 135             wi = 0.0;
 136
 137             for (i3 = 1; i3 <= ifp1; i3 += ip1)
 138             {
 139                 for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
 140                 {
 141                     for (i2 = i1; i2 <= ip3; i2 += ifp2)
 142                     {
 143                         k1 = i2;
 144                         k2 = k1 + ifp1;
 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;
 149                         data[k1] += tempr;
 150                         data[k1 + 1] += tempi;
 151                     }
 152                 }
 153
 154                 wr = (wtemp = wr)*wpr - wi*wpi + wr;
 155                 wi = wi*wpr + wtemp*wpi + wi;
 156             }
 157             ifp1 = ifp2;
 158         }
 159         nprev *= n;
 160     }
 161
 162
 163     // if forward transform : renumber after transform
 164
 165     if (isign == FORWARD_TRANSFORM)
 166     {
 167         fftRenumber(field, nn);
 168     }
 169
 170
 171     // scale result (symmetric scale both forward and inverse transform)
 172
 173     scalar recRootN = 1.0/sqrt(scalar(ntot));
 174
 175     forAll(field, i)
 176     {
 177         field[i] *= recRootN;
 178     }
 179 }
 180
 181
 182 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 183
 184 #undef SWAP
 185 #undef TWOPI
 186
 187 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 188
 189 tmp<complexField> fft::forwardTransform
 190 (
 191     const tmp<complexField>& tfield,
 192     const labelList& nn
 193 )
 194 {
 195     tmp<complexField> tfftField(new complexField(tfield));
 196
 197     transform(tfftField(), nn, FORWARD_TRANSFORM);
 198
 199     tfield.clear();
 200
 201     return tfftField;
 202 }
 203
 204
 205 tmp<complexField> fft::reverseTransform
 206 (
 207     const tmp<complexField>& tfield,
 208     const labelList& nn
 209 )
 210 {
 211     tmp<complexField> tifftField(new complexField(tfield));
 212
 213     transform(tifftField(), nn, REVERSE_TRANSFORM);
 214
 215     tfield.clear();
 216
 217     return tifftField;
 218 }
 219
 220
 221 tmp<complexVectorField> fft::forwardTransform
 222 (
 223     const tmp<complexVectorField>& tfield,
 224     const labelList& nn
 225 )
 226 {
 227     tmp<complexVectorField> tfftVectorField
 228     (
 229         new complexVectorField
 230         (
 231             tfield().size()
 232         )
 233     );
 234
 235     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 236     {
 237         tfftVectorField().replace
 238         (
 239             cmpt,
 240             forwardTransform(tfield().component(cmpt), nn)
 241         );
 242     }
 243
 244     tfield.clear();
 245
 246     return tfftVectorField;
 247 }
 248
 249
 250 tmp<complexVectorField> fft::reverseTransform
 251 (
 252     const tmp<complexVectorField>& tfield,
 253     const labelList& nn
 254 )
 255 {
 256     tmp<complexVectorField> tifftVectorField
 257     (
 258         new complexVectorField
 259         (
 260             tfield().size()
 261         )
 262     );
 263
 264     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 265     {
 266         tifftVectorField().replace
 267         (
 268             cmpt,
 269             reverseTransform(tfield().component(cmpt), nn)
 270         );
 271     }
 272
 273     tfield.clear();
 274
 275     return tifftVectorField;
 276 }
 277
 278
 279 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 280
 281 } // End namespace Foam
 282
 283 // ************************************************************************* //