Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/map_2d.hh
e96d4bebeabcfeatures/netcdf-mpi

Last Change
Raw File

map_2d.hh
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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __MAP_2D_HH__
#define __MAP_2D_HH__
/* -------------------------------------------------------------------------- */
#include <vector>
#include <iostream>
#include <fftw3.h>
#include "tamaas.hh"
#include <fstream>
#include "my_vector.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
template <typename T>
class Map2d {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
// class iterator;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Map2d(UInt n1, UInt n2, Real L0 = 1., Real L1 = 1.);
Map2d(const UInt * n, const Real * L);
virtual ~Map2d();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <bool periodic>
std::vector<std::pair<UInt,UInt> > getNeighborIndexes(UInt i, UInt j);
template <bool periodic>
std::vector<std::pair<UInt,UInt> > getNeighborIndexes(std::pair<UInt,UInt> i);
// //! return an iterator over the neighbors
// iterator beginNeighborhood();
// //! return the end iterator over the neighbors
// iterator endNeighborhood();
//! generate fft plan
void generateFFTPlan();
//! destroy fft plan
void destroyFFTPlan();
//! perform fft transform
void FFTTransform(Map2d<std::complex<Real> > & output,int nthreads = 1);
//! perform inverse fft transform
void FFTITransform(Map2d<std::complex<Real> > & output, int nthreads = 1);
//! perform inverse fft transform
void FFTITransform(Map2d<Real> & output, int nthreads = 1);
//! perform fft transform
void FFTTransform(int nthreads = 1);
// //! inverse fft transform
void FFTITransform(int nthreads = 1);
//compute compute component by component product and put result in object
template <typename T1, typename T2>
void fullProduct(Map2d<T1> & s1,Map2d<T2> &s2);
//! set the grid size
void setGridSize(const UInt n[2]);
//! set the grid size
void setGridSize(UInt n0,UInt n1);
//! rescale height profile
void rescaleHeights(Real factor);
//! Write heights array to text file usable with plot programs
void dumpToTextFile(std::string filename);
// set all the members of the map to a given value
void operator=(const T & e);
template <template<typename> class surf, typename T2>
void operator=(const surf<T2> & s);
template <template<typename> class surf, typename T2>
void operator=(surf<T2> & s);
/* -------------------------------------------------------------------------- */
private:
//! write an entry to a stringstream
void writeValue(std::ostream & os, const T & val);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
// access to the i,j-th value in the surface
inline void operator+=(const Map2d<T> & s);
// access to the i,j-th value in the surface
inline void operator+=(const T & val);
// access to the i,j-th value in the surface
inline void operator-=(const T & val);
// access to the i,j-th value in the surface
inline T & operator()(UInt i,UInt j);
// access to the i,j-th value in the surface
inline T & operator()(std::pair<UInt,UInt> i);
// access to the i,j-th value in the surface
inline const T & operator()(UInt i,UInt j)const;
// access to the i,j-th value in the surface
inline const T & operator()(std::pair<UInt,UInt> i)const;
// access to the i,j-th value in the surface
inline T & at(UInt i,UInt j);
// access to the i-th value in the array
inline T & at(UInt i);
// access to the i,j-th value in the surface
inline T & at(std::pair<UInt,UInt> i);
// access to the i,j-th value in the surface
inline const T & at(UInt i,UInt j) const;
// access to the i,j-th value in the surface
inline const T & at(std::pair<UInt,UInt> i) const;
// access to the i,j-th value in the surface
inline const T & at(UInt i) const;
// access to the i-th value in the array
inline T & operator()(int i);
// access to the i-th value in the array
// inline T & operator[](std::pair<UInt,UInt> i);
// access to the i-th value in the array
inline const T & operator()(int i) const;
//! get lateral number of discretization points
inline UInt size(UInt dir)const {return n[dir];};
//! get both lateral number of discretization points
inline const UInt* sizes()const {return n;};
//! get the lateral size of the sample
inline const Real * getLs() const {return L;};
//! get the lateral size of the sample
inline Real getL(UInt dir) const {return L[dir];};
//! get the lateral size of the sample
inline void setL(Real L[2]){this->L[0] = L[0];this->L[1] = L[1];};
//! get the lateral size of the sample
inline void setL(Real L0,Real L1){this->L[0] = L0;this->L[1] = L1;};
const T * getInternalData() const{return data.getPtr();};
//! Write heights array to text file usable with plot programs
template <bool consider_zeros=true>
void dumpToTextFile(std::string filename) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! number of points in each direction
UInt n[2];
//! size in each direction
Real L[2];
//! internal storage of data
MyVector<T> data;
//! fftw plan
fftw_plan plan;
};
/* -------------------------------------------------------------------------- */
//#define _at(I,J,N) I+J*N[1]
#define _at(I,J,N) I*N[1]+J
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::at(UInt i,UInt j){
if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
<< i << " " << j << " "
<< n[0] << " " << n[1] << " "
<< _at(i,j,n));
return data[_at(i,j,n)];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::at(UInt i,UInt j)const {
if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
<< i << " " << j << " "
<< n[0] << " " << n[1] << " "
<< _at(i,j,n));
return data[_at(i,j,n)];
}
/* -------------------------------------------------------------------------- */
#undef _at
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::at(UInt i){
return data[i];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::at(UInt i) const{
return data[i];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::operator()(UInt i,UInt j){
return this->at(i,j);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::at(std::pair<UInt,UInt> i){
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::at(std::pair<UInt,UInt> i)const {
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::operator()(std::pair<UInt,UInt> i){
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::operator()(UInt i,UInt j)const{
return this->at(i,j);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::operator()(std::pair<UInt,UInt> i)const {
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::operator()(int i){
return data[i];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::operator()(int i) const{
return data[i];
}
/* -------------------------------------------------------------------------- */
// template <typename T>
// inline T & Map2d<T>::operator[](std::pair<UInt,UInt> c){
// return this->at(c.first,c.second);
// }
// /* -------------------------------------------------------------------------- */
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>
void Map2d<T>::operator+=(const T & e){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) += e;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator+=(const Map2d<T> & s){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) += s.at(i);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <bool periodic>
std::vector<std::pair<UInt,UInt> > Map2d<T>::getNeighborIndexes(UInt i, UInt j){
std::vector<std::pair<UInt,UInt> > vres;
std::pair<int,int> res;
for (int ii = -1; ii < 2; ++ii) {
res.first = ii+static_cast<int>(i);
if (res.first < 0 && !periodic) continue;
if (res.first >= (int)n[0] && !periodic) continue;
if (res.first < 0 && periodic) res.first += n[0];
if (res.first >= (int)n[0] && periodic) res.first -= n[0];
for (int jj = -1; jj < 2; ++jj) {
res.second = jj+static_cast<int>(j);
if (res.second < 0 && !periodic) continue;
if (res.second >= (int)n[1] && !periodic) continue;
if (res.second < 0 && periodic) res.second += n[1];
if (res.second >= (int)n[1] && periodic) res.second -= n[1];
if (ii == 0 && jj == 0) continue;
{
UInt i = res.first;
UInt j = res.second;
if (i+j*n[0] >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
<< i << " " << j << " "
<< n[0] << " " << n[1] << " "
<< i+j*n[0]);
}
vres.push_back(res);
}
}
return vres;
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <bool periodic>
std::vector<std::pair<UInt,UInt> > Map2d<T>
::getNeighborIndexes(std::pair<UInt,UInt> i){
return this->getNeighborIndexes<periodic>(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <typename T1, typename T2>
void Map2d<T>::fullProduct(Map2d<T1> & s1,Map2d<T2> & s2){
auto n1 = s1.sizes();
auto n2 = s2.sizes();
if (n1[0] != n2[0] ||
n1[1] != n2[1]) SURFACE_FATAL("fullproduct of non compatible matrices");
if (n[0] != n1[0] ||
n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices")
for (UInt i = 0 ; i < n[0] ; ++i)
for (UInt j = 0 ; j < n[1] ; ++j){
this->at(i,j) = s1(i,j)* s2(i,j);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <bool consider_zeros>
void Map2d<T>::dumpToTextFile(std::string filename) const {
UInt m = this->n[0];
UInt n = this->n[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 < m; i++) {
for (UInt j=0; j< n; j++) {
if (!consider_zeros &&
std::norm(this->at(i,j)) < zero_threshold*zero_threshold)
continue;
file << i << " " << j << " ";
file << std::scientific << this->at(i,j);
file << std::endl;
}
file << std::endl;
}
file.close();
}
// template <typename T>
// class Map2d<T>::iterator {
// public:
// iterator(std::vector<T> & data,
// std::queue<std::pair<UInt,UInt> > & indices):
// data(data),indices(indices){
// }
// inline void operator++(){
// indices.pop();
// }
// inline T & operator *(){
// return data[indices.front];
// }
// bool operator != (Map2d<T>::iterator & it){
// return (indices.front == it.indices.front);
// }
// private:
// std::vector<T> & data;
// std::queue<std::pair<UInt,UInt> > indices;
// bool finished;
// };
/* -------------------------------------------------------------------------- */
template <typename T>
template <template<typename> class surf, typename T2>
void Map2d<T>::operator=(const surf<T2> & s){
this->setGridSize(s.sizes());
UInt m = this->sizes()[0];
UInt n = this->sizes()[1];
#pragma omp parallel for
for (UInt i = 0; i < m*n; i++) {
this->at(i) = s(i);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <template<typename> class surf, typename T2>
void Map2d<T>::operator=(surf<T2> & s){
this->setGridSize(s.sizes());
UInt m = this->sizes()[0];
UInt n = this->sizes()[1];
#pragma omp parallel for
for (UInt i = 0; i < m*n; i++) {
this->at(i) = s(i);
}
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif /* __MAP_2D_HH__ */
c4science · Help