Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NekFFTW.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NekFFTW.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Wrapper around FFTW library
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
39 
40 namespace Nektar
41 {
42  namespace LibUtilities
43  {
44  string NekFFTW::className
47 
49  : NektarFFT(N)
50  {
51  m_wsp = Array<OneD, NekDouble>(m_N);
52  phys = Array<OneD,NekDouble>(m_N);
53  coef = Array<OneD,NekDouble>(m_N);
54 
55  plan_forward = fftw_plan_r2r_1d(m_N, &phys[0], &coef[0],
56  FFTW_R2HC, FFTW_ESTIMATE);
57  plan_backward = fftw_plan_r2r_1d(m_N, &coef[0], &phys[0],
58  FFTW_HC2R, FFTW_ESTIMATE);
59 
60  m_FFTW_w = Array<OneD,NekDouble>(m_N);
61  m_FFTW_w_inv = Array<OneD,NekDouble>(m_N);
62 
63  m_FFTW_w[0] = 1.0/(NekDouble)m_N;
64  m_FFTW_w[1] = 0.0;
65 
66  m_FFTW_w_inv[0] = m_N;
67  m_FFTW_w_inv[1] = 0.0;
68 
69  for(int i=2;i<m_N;i++)
70  {
71  m_FFTW_w[i] = m_FFTW_w[0];
72  m_FFTW_w_inv[i] = m_FFTW_w_inv[0];
73  }
74  }
75 
76  // Distructor
78  {
79 
80  }
81 
82  // Forward transformation
84  Array<OneD,NekDouble> &inarray,
85  Array<OneD,NekDouble> &outarray)
86  {
87  Vmath::Vcopy(m_N, inarray, 1, phys, 1);
88 
89  fftw_execute(plan_forward);
90 
92 
93  Vmath::Vcopy(m_N, coef, 1, outarray, 1);
94  }
95 
96  // Backward transformation
98  Array<OneD,NekDouble> &inarray,
99  Array<OneD,NekDouble> &outarray)
100  {
101  Vmath::Vcopy(m_N, inarray, 1, coef, 1);
102 
104 
105  fftw_execute(plan_backward);
106 
107  Vmath::Vcopy(m_N, phys, 1, outarray, 1);
108  }
109 
110  // Reshuffle FFTW2Nek
111  void NekFFTW::Reshuffle_FFTW2Nek(Array<OneD,NekDouble> &coef)
112  {
113  int halfN = m_N/2;
114 
115  m_wsp[1] = coef[halfN];
116 
117  Vmath::Vcopy(halfN, coef, 1, m_wsp, 2);
118 
119  for(int i = 0; i < (halfN - 1); i++)
120  {
121  m_wsp[(m_N-1)-2*i] = coef[halfN+1+i];
122  }
123 
124  Vmath::Vmul(m_N, m_wsp, 1, m_FFTW_w, 1, coef, 1);
125 
126  return;
127  }
128 
129  // Reshuffle Nek2FFTW
130  void NekFFTW::Reshuffle_Nek2FFTW(Array<OneD,NekDouble> &coef)
131  {
132  int halfN = m_N/2;
133 
134  Vmath::Vmul(m_N, coef, 1, m_FFTW_w_inv, 1, coef, 1);
135 
136  m_wsp[halfN] = coef[1];
137 
138  Vmath::Vcopy(halfN, coef, 2, m_wsp, 1);
139 
140  for(int i = 0; i < (halfN-1); i++)
141  {
142  m_wsp[halfN+1+i] = coef[(m_N-1)-2*i];
143  }
144 
145  Vmath::Smul(m_N, 1.0/m_N, m_wsp, 1, coef, 1);
146 
147  return;
148  }
149  }
150 }