Homec4science
Diffusion tamaas (master) Browse

tamaas/src/fftransform.cpp
9251b566b1bfmaster

Last Change
Raw File

fftransform.cpp
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "fftransform.hh"
__BEGIN_TAMAAS__
template<typename T>
FFTransform<T>::FFTransform(UInt size,
Surface<T> & real,
SurfaceComplex<T> & spectral):
real(real),
spectral(spectral),
forward_plan(NULL),
backward_plan(NULL),
size(size),
reduced_size(size / 2 + 1)
{
FFTFlag flag = ESTIMATE;
spectral_buffer = new fftw_complex[size * reduced_size];
Real * in = const_cast<Real *>(real.getInternalData());
forward_plan = fftw_plan_dft_r2c_2d(size, size,
in,
spectral_buffer,
flag);
backward_plan = fftw_plan_dft_c2r_2d(size, size,
spectral_buffer,
in,
flag);
}
/* -------------------------------------------------------------------------- */
template<typename T>
FFTransform<T>::~FFTransform() {
fftw_destroy_plan(forward_plan);
fftw_destroy_plan(backward_plan);
delete[] spectral_buffer;
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransform<T>::forwardTransform() {
fftw_execute(forward_plan);
writeBuffer();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransform<T>::backwardTransform() {
loadBuffer();
fftw_execute(backward_plan);
normalize();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransform<T>::normalize() {
#pragma omp parallel for
for (UInt i = 0 ; i < size * size ; ++i)
real(i) /= size * size;
}
/* -------------------------------------------------------------------------- */
template <typename T>
void FFTransform<T>::writeBuffer() {
#pragma omp parallel
{
#pragma omp for
// First line of spectral surface is different
for (UInt j = 1 ; j < reduced_size-1 ; j++) {
std::complex<T> z(spectral_buffer[j][0],
spectral_buffer[j][1]);
spectral(0, j) = z;
spectral(0, size-j) = std::conj(z);
}
#pragma omp for
// First column just needs copying
for (UInt i = 0 ; i < size ; ++i) {
UInt index = i * reduced_size;
std::complex<T> z(spectral_buffer[index][0],
spectral_buffer[index][1]);
spectral(i, 0) = z;
}
#pragma omp for
// Copy reduced_size-1 column
for (UInt i = 0 ; i < size ; ++i) {
UInt index = i * reduced_size + reduced_size-1;
std::complex<T> z(spectral_buffer[index][0],
spectral_buffer[index][1]);
spectral(i, reduced_size-1) = z;
}
#pragma omp for collapse(2)
// Rest of spectral surface is entertwined
for (UInt i = 1 ; i < size ; i++) {
for (UInt j = 1 ; j < reduced_size-1 ; j++) {
UInt index = i * reduced_size + j;
std::complex<T> z(spectral_buffer[index][0],
spectral_buffer[index][1]);
spectral(i, j) = z;
spectral(size-i, size-j) = std::conj(z);
}
}
} // end #pramga omp parallel
}
/* -------------------------------------------------------------------------- */
template <typename T>
void FFTransform<T>::loadBuffer() {
#pragma omp parallel for collapse(2)
for (UInt i = 0 ; i < size ; i++) {
for (UInt j = 0 ; j < reduced_size ; j++) {
UInt index = i * reduced_size + j;
spectral_buffer[index][0] = spectral(i, j).real();
spectral_buffer[index][1] = spectral(i, j).imag();
}
}
}
/* -------------------------------------------------------------------------- */
template class FFTransform<Real>;
__END_TAMAAS__
c4science · Help