Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/fftransform_fftw.cpp
51a6f524e10afeatures/declarative

Last Change
Raw File

fftransform_fftw.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_fftw.hh"
__BEGIN_TAMAAS__
template<typename T>
FFTransformFFTW<T>::FFTransformFFTW(UInt size,
Surface<T> & real,
SurfaceComplex<T> & spectral):
FFTransform<T>(size, real, spectral),
forward_plan(NULL),
backward_plan(NULL)
{
this->spectral_buffer = new fftw_complex[this->size * this->reduced_size];
Real * in = const_cast<Real *>(this->real.getInternalData());
forward_plan = fftw_plan_dft_r2c_2d(this->size, this->size,
in,
this->spectral_buffer,
FFTW_ESTIMATE);
backward_plan = fftw_plan_dft_c2r_2d(this->size, this->size,
this->spectral_buffer,
in,
FFTW_ESTIMATE);
}
/* -------------------------------------------------------------------------- */
template<typename T>
FFTransformFFTW<T>::~FFTransformFFTW() {
fftw_destroy_plan(forward_plan);
fftw_destroy_plan(backward_plan);
delete[] this->spectral_buffer;
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransformFFTW<T>::forwardTransform() {
fftw_execute(forward_plan);
writeBuffer();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransformFFTW<T>::backwardTransform() {
loadBuffer();
fftw_execute(backward_plan);
this->normalize();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void FFTransformFFTW<T>::writeBuffer() {
#pragma omp parallel
{
#pragma omp for
// First line of this->spectral surface is different
for (UInt j = 1 ; j < this->reduced_size-1 ; j++) {
std::complex<T> z(this->spectral_buffer[j][0],
this->spectral_buffer[j][1]);
this->spectral(0, j) = z;
this->spectral(0, this->size-j) = std::conj(z);
}
#pragma omp for
// First column just needs copying
for (UInt i = 0 ; i < this->size ; ++i) {
UInt index = i * this->reduced_size;
std::complex<T> z(this->spectral_buffer[index][0],
this->spectral_buffer[index][1]);
this->spectral(i, 0) = z;
}
#pragma omp for
// Copy reduced_size-1 column
for (UInt i = 0 ; i < this->size ; ++i) {
UInt index = i * this->reduced_size + this->reduced_size-1;
std::complex<T> z(this->spectral_buffer[index][0],
this->spectral_buffer[index][1]);
this->spectral(i, this->reduced_size-1) = z;
}
#pragma omp for collapse(2)
// Rest of this->spectral surface is entertwined
for (UInt i = 1 ; i < this->size ; i++) {
for (UInt j = 1 ; j < this->reduced_size-1 ; j++) {
UInt index = i * this->reduced_size + j;
std::complex<T> z(this->spectral_buffer[index][0],
this->spectral_buffer[index][1]);
this->spectral(i, j) = z;
this->spectral(this->size-i, this->size-j) = std::conj(z);
}
}
} // end #pramga omp parallel
}
/* -------------------------------------------------------------------------- */
template <typename T>
void FFTransformFFTW<T>::loadBuffer() {
#pragma omp parallel for collapse(2)
for (UInt i = 0 ; i < this->size ; i++) {
for (UInt j = 0 ; j < this->reduced_size ; j++) {
UInt index = i * this->reduced_size + j;
this->spectral_buffer[index][0] = this->spectral(i, j).real();
this->spectral_buffer[index][1] = this->spectral(i, j).imag();
}
}
}
/* -------------------------------------------------------------------------- */
template class FFTransformFFTW<Real>;
__END_TAMAAS__
c4science · Help