src/randomProcesses/fft/fftRenumber.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 Description
  25     Multi-dimensional renumbering used in the Numerical Recipes
  26    fft routine. This version is recursive, so works in n-d :
  27    determined by the length of array nn
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "fftRenumber.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 // recursively evaluate the indexing necessary to do the folding of
  41 // the fft data. We recurse until we have the indexing ncessary for
  42 // the folding in all directions.
  43
  44 void fftRenumberRecurse
  45 (
  46     List<complex>& data,
  47     List<complex>& renumData,
  48     const labelList& nn,
  49     label nnprod,
  50     label ii,
  51     label l1,
  52     label l2
  53 )
  54 {
  55     if (ii == nn.size())
  56     {
  57         // we've worked out the renumbering scheme. Now copy
  58         // the components across
  59
  60         data[l1] = complex(renumData[l2].Re(),renumData[l2].Im());
  61     }
  62     else
  63     {
  64         // do another level of folding. First work out the
  65         // multiplicative value of the index
  66
  67         nnprod /= nn[ii];
  68         label i_1(0);
  69
  70         for (label i=0; i<nn[ii]; i++)
  71         {
  72             // now evaluate the indices (both from array 1 and to
  73             // array 2). These get multiplied by nnprod to (cumulatively)
  74             // find the real position in the list corresponding to
  75             // this set of indices.
  76
  77             if (i<nn[ii]/2)
  78             {
  79                 i_1 = i + nn[ii]/2;
  80             }
  81             else
  82             {
  83                 i_1 = i - nn[ii]/2;
  84             }
  85
  86
  87             // go to the next level of recursion.
  88
  89             fftRenumberRecurse
  90             (
  91                 data,
  92                 renumData,
  93                 nn,
  94                 nnprod,
  95                 ii+1,
  96                 l1+i*nnprod,
  97                 l2+i_1*nnprod
  98             );
  99         }
 100     }
 101 }
 102
 103
 104 // fftRenumber : fold the n-d data array to get the fft components in
 105 // the right places.
 106
 107 #include "fftRenumber.H"
 108
 109 void fftRenumber
 110 (
 111     List<complex>& data,
 112     const labelList& nn
 113 )
 114 {
 115     List<complex> renumData(data);
 116
 117     label nnprod(1);
 118     for (label i=0; i<nn.size(); i++)
 119     {
 120         nnprod *= nn[i];
 121     }
 122
 123     label ii(0), l1(0), l2(0);
 124
 125     fftRenumberRecurse
 126     (
 127         data,
 128         renumData,
 129         nn,
 130         nnprod,
 131         ii,
 132         l1,
 133         l2
 134     );
 135 }
 136
 137
 138 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 139
 140 } // End namespace Foam
 141
 142 // ************************************************************************* //