Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/map_2d.cpp
7d473b01cbbamaxwell_viscoelastic…

Last Change
Raw File

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 <fstream>
#include <sstream>
#include <fftw3.h>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::rescaleHeights(Real C){
/* compute actual RMS */
for (UInt i = 0 ; i < n[0] ; ++i)
for (UInt j = 0 ; j < n[1] ; ++j){
this->at(i,j) *= C;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::setGridSize(UInt n0, UInt n1){
n[0] = n0;
n[1] = n1;
data.resize(n[0]*n[1]);
data.assign(n[0]*n[1],T(0.));
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::setGridSize(const UInt n[2]){
setGridSize(n[0],n[1]);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator=(const T & e){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) = 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;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::Map2d(const UInt * n, const Real *L) {
setGridSize(n[0],n[1]);
this->L[0] = L[0];
this->L[1] = L[1];
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::~Map2d() {
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTTransform(int nthreads){
SURFACE_FATAL("impossible situation");
throw;
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTTransform(Map2d<std::complex<Real> > & output, int nthreads){
output = *this;
output.FFTTransform(nthreads);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(Map2d<std::complex<Real> > & output, int nthreads){
output = *this;
output.FFTITransform(nthreads);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(Map2d<Real> & output, int nthreads){
const UInt * n = this->sizes();
const Real * Ls = this->getLs();
Map2d<std::complex<Real> > _output(n,Ls);
_output = *this;
_output.FFTITransform(nthreads);
for (UInt i = 0; i < n[0]*n[1]; i++) {
output(i) = _output(i).real();
}
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::FFTTransform(int nthreads){
STARTTIMER("initFFTW");
fftw_init_threads();
fftw_plan_with_nthreads(nthreads);
fftw_complex * in = (fftw_complex*)&data[0];
fftw_plan p = fftw_plan_dft_2d(n[0],n[1], in,in,FFTW_FORWARD, FFTW_ESTIMATE);
STOPTIMER("initFFTW");
STARTTIMER("FFTW");
fftw_execute(p);
STOPTIMER("FFTW");
STARTTIMER("closeFFTW");
fftw_destroy_plan(p);
fftw_cleanup_threads();
STOPTIMER("closeFFTW");
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(int nthreads){
SURFACE_FATAL("impossible situation");
throw;
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::FFTITransform(int nthreads){
for (UInt i = 0; i < n[0]*n[1]; ++i) {
data[i] = std::conj(data[i]);
}
FFTTransform(nthreads);
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>;
c4science · Help