Page MenuHomec4science
Files F67578099

map_2d.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Jun 23, 04:36

map_2d.cpp
View Options

/**
*
* @author Guillaume Anciaux <guillaume.anciaux@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 "map_2d.hh"
#include "fft_plan_manager.hh"
#include <fstream>
#include <sstream>
#include <fftw3.h>
#include <omp.h>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
template <typename T>
void Map2d<T>::rescaleHeights(Real C){
/* compute actual RMS */
for (UInt i = 0 ; i < this->n[0] ; ++i)
for (UInt j = 0 ; j < this->n[1] ; ++j){
this->at(i,j) *= C;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::setGridSize(UInt n0, UInt n1){
const UInt n[] = {n0, n1};
this->resize(n);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::setGridSize(const UInt n[2]){
this->resize(n);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator=(const T & e){
this->Grid<T, 2>::operator=(e);
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::Map2d(UInt n0, UInt n1, Real L0, Real L1) {
setGridSize(n0,n1);
this->L[0] = L0;
this->L[1] = L1;
this->plan = NULL;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::Map2d(const UInt * n, const Real *L) {
setGridSize(n);
this->L[0] = L[0];
this->L[1] = L[1];
this->plan = NULL;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::~Map2d() {
destroyFFTPlan();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::generateFFTPlan() {
// this should not do anything
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::generateFFTPlan() {
fftw_complex * in = (fftw_complex*) &data[0];
this->plan = fftw_plan_dft_2d(n[0], n[1], in, in, FFTW_FORWARD, FFTW_ESTIMATE);
}
/* -------------------------------------------------------------------------- */
template<typename T>
void Map2d<T>::destroyFFTPlan() {
if (this->plan) {
fftw_destroy_plan(this->plan);
this->plan = NULL;
}
FFTPlanManager::get().destroyPlan(this->getInternalData());
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTTransform(){
SURFACE_FATAL("impossible situation");
throw;
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTTransform(Map2d<std::complex<Real> > & output)const{
output = *this;
output.FFTTransform();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(Map2d<std::complex<Real> > & output)const{
output = *this;
output.FFTITransform();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(Map2d<Real> & output) const{
const UInt * n = this->sizes();
const Real * Ls = this->getLs();
Map2d<std::complex<Real> > _output(n,Ls);
_output = *this;
_output.FFTITransform();
for (UInt i = 0; i < n[0]*n[1]; i++) {
output(i) = _output(i).real();
}
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::FFTTransform(){
STARTTIMER("FFTW");
if (!this->plan)
generateFFTPlan();
fftw_execute(this->plan);
STOPTIMER("FFTW");
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(){
SURFACE_FATAL("impossible situation");
throw;
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::FFTITransform(){
#pragma omp parallel for
for (UInt i = 0; i < n[0]*n[1]; ++i) {
data[i] = std::conj(data[i]);
}
FFTTransform();
#pragma omp parallel for
for (UInt i = 0; i < n[0]*n[1]; ++i) {
data[i] = std::conj(data[i]);
data[i] *= 1./(n[0]*n[1]);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::dumpToTextFile(std::string filename) {
UInt n1 = this->size(0);
UInt n2 = this->size(1);
std::cout << "Writing surface file " << filename << std::endl;
std::ofstream file(filename.c_str());
if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename);
file.precision(15);
for (UInt i = 0; i < n1; i++) {
for (UInt j=0; j< n2; j++) {
file << i << " " << j << " ";
this->writeValue(file,this->at(i,j));
file << std::endl;
}
file << std::endl;
}
file.close();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::writeValue(std::ostream & os, const T & val) {
os << val;
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::writeValue(std::ostream & os, const complex & val) {
os << std::scientific << val.real();
os << " " << val.imag();
}
/* -------------------------------------------------------------------------- */
template class Map2d<complex>;
template class Map2d<Real>;
template class Map2d<int>;
template class Map2d<unsigned int>;
template class Map2d<bool>;
template class Map2d<unsigned long>;
__END_TAMAAS__

Event Timeline

c4science · Help