src/randomProcesses/fft/fft.C

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