Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/core/surface.cpp
cb641f5bb816features/netcdf-mpi

Last Change
Raw File

surface.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include <fstream>
#include <sstream>
#ifdef USE_CUDA
#include <thrust/random.h>
#else
#include <random>
#endif
/* -------------------------------------------------------------------------- */
#include "surface.hh"
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
#ifdef USING_IOHELPER
#include "dumper_paraview.hh"
#endif
/* -------------------------------------------------------------------------- */
namespace tamaas {
std::string fname = "surface_height.pdf";
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::expandByLinearInterpolation() {
UInt n = this->size();
auto old_s = std::make_unique<Surface>(n, 1.);
int new_n = n * 2;
*old_s = *this;
this->setGridSize(new_n);
for (UInt i = 0; i < n; i++) {
for (UInt j = 0; j < n; j++) {
(*this)(i, j) = (*old_s)(i, j);
}
}
}
Real sinc(Real x) { return sin(M_PI * x) / (M_PI * x); }
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::expandByShannonInterpolation(int /*factor*/) {
/* /!\ this->FFTTransform(); /!\ */
TAMAAS_EXCEPTION("Shannon interpolation needs old FFTTransform interface");
#if 0
UInt n = this->size();
UInt N = std::pow(2,factor)*this->size();
std::cout << N << " " << factor << std::endl;
Surface s = Surface(N,1.);
for (UInt i = 0; i < n/2; ++i)
for (UInt j = 0; j < n/2; ++j) {
s(i,j) = (*this)(i,j);
}
for (UInt i = n/2; i < n; ++i)
for (UInt j = 0; j < n/2; ++j) {
s(N-n+i,j) = (*this)(i,j);
}
for (UInt i = n/2; i < n; ++i)
for (UInt j = n/2; j < n; ++j) {
s(N-n+i,N-n+j) = (*this)(i,j);
}
for (UInt i = 0; i < n/2; ++i)
for (UInt j = n/2; j < n; ++j) {
s(i,N-n+j) = (*this)(i,j);
}
// s.dumpToTextFile<false,true>("temp.csv");
/* /!\ s.FFTITransform(); /!\ */
this->setGridSize(N);
Real f = std::pow(2,factor);
f = f*f;
for (UInt i = 0; i < N*N; ++i)
(*this)(i) = f*s(i);
// int new_n = 2*n;
// Surface s(*this);
// this->setGridSize(new_n);
// // put the original signal with zeros for to be interpolated points
// for (int i=0; i<n; i++)
// for (int j=0; j<n; j++)
// (*this)(2*i,j*2) = s(i,j);
// // do shannon interpolation
// int rcut = 64;
// //interpolate lines
// for (int i=0; i<new_n; i++) {
// for (int j=0; j<new_n; j++) {
// if (i%2 == 0) continue;
// // Real loop on local neighbors to interpolate
// for (int k=-rcut; k<rcut+1; k++) {
// int i2 = i + k;
// if (i2 < 0) i2 += new_n;
// if (i2 >= new_n) i2 -= new_n;
// if (i2%2) continue;
// Real d = fabs(k);
// s[i+j*new_n] += s(i2,j)*sinc(d/2);
// }
// }
// }
// //interpolate columns
// for (int i=0; i<new_n; i++) {
// for (int j=0; j<new_n; j++) {
// if (j%2 == 0) continue;
// // Real loop on local neighbors to interpolate
// for (int k=-rcut; k<rcut+1; k++) {
// int j2 = j + k;
// if (j2 < 0) j2 += new_n;
// if (j2 >= new_n) j2 -= new_n;
// if (j2%2) continue;
// Real d = fabs(k);
// s[i+j*new_n] += complex(s[i+j2*new_n].real()*sinc(d/2),
// s[i+j2*new_n].imag()*sinc(d/2));
// }
// }
// }
// // for (int i=0; i<n; i++) {
// // for (int j=0; j<n; j++) {
// // s[2*i+j*2*new_n].real() = 0;
// // s[2*i+j*2*new_n].im = 0;
// // }
// // }
#endif
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::expandByMirroring() {
UInt n = this->size();
UInt new_n = n * 2;
Surface s(*this);
this->setGridSize(new_n);
for (UInt i = 0; i < new_n; i++) {
for (UInt j = 0; j < new_n; j++) {
UInt i1 = i;
if (i1 > n)
i1 = new_n - i1;
UInt j1 = j;
if (j1 > n)
j1 = new_n - j1;
if (i == n || j == n) {
(*this)(i, j) = 0.;
continue;
}
(*this)(i, j) = s(i1, j1);
}
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::expandByPuttingZeros() {
UInt n = this->size();
UInt new_n = n * 2;
Surface s(*this);
this->setGridSize(new_n);
for (UInt i = 0; i < n; i++) {
for (UInt j = 0; j < n; j++) {
(*this)(2 * i, j * 2) = s(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::contractByTakingSubRegion() {
UInt n = this->size();
Surface s(*this);
this->setGridSize(n / 2);
for (UInt i = 0; i < n / 2; i++) {
for (UInt j = 0; j < n / 2; j++) {
(*this)(i, j) = s(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::contractByIgnoringMidPoints() {
UInt n = this->size();
Surface s(*this);
this->setGridSize(n / 2);
for (UInt i = 0; i < n; i++) {
for (UInt j = 0; j < n; j++) {
if (i % 2)
continue;
if (j % 2)
continue;
(*this)(i / 2, j / 2) = s(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
// void Surface<T>::dumpToTextFile(string filename, bool pr) {
// cout << "Writing surface file " << filename << endl;
// ofstream file(filename.c_str());
// file.precision(15);
// if (pr)
// {
// for (int i=0; i<n; i++) {
// for (int j=0; j<n; j++) {
// file << i << " " << j << " ";
// if (heights[i+j*n].real() < 1e-10)
// file << scientific << -1e10 << " " << 0;
// else
// file << scientific << heights[i+j*n].real() << " " <<
// heights[i+j*n].im;
// file << endl;
// }
// file << endl;
// }
// }
// else
// {
// for (int i=0; i<n; i++) {
// for (int j=0; j<n; j++) {
// file << i << " " << j << " ";
// file << scientific << heights[i+j*n].real() << " " <<
// heights[i+j*n].im; file << endl;
// }
// }
// }
// file.close();
// }
/* -------------------------------------------------------------------------- */
// Surface * Surface<T>::loadFromTextFile(std::string filename) {
// std::string line;
// std::ifstream myfile (filename.c_str());
// if (!myfile.is_open()){
// SURFACE_FATAL("cannot open file " << filename);
// }
// Surface * my_surface = nullptr;
// UInt n = 0;
// if (myfile.is_open())
// {
// while (! myfile.eof() )
// {
// getline (myfile,line);
// if (line[0] == '#')
// continue;
// if (line.size() == 0)
// continue;
// std::stringstream s;
// s << line;
// int i,j;
// s >> i;
// s >> j;
// if (n < i) n = i;
// if (n < j) n = j;
// }
// ++n;
// std::cout << "read " << n << " x " << n << " grid points" << std::endl;
// myfile.clear(); // forget we hit the end of file
// myfile.seekg(0, std::ios::beg); // move to the start of the file
// my_surface = new Surface(n);
// while (! myfile.eof() )
// {
// getline (myfile,line);
// if (line[0] == '#')
// continue;
// if (line.size() == 0)
// continue;
// std::stringstream s;
// s << line;
// int i,j;
// Real h,hIm;
// s >> i;
// s >> j;
// s >> h;
// s >> hIm;
// my_surface->heights[i+j*n] = complex(h,hIm);
// }
// myfile.close();
// }
// return my_surface;
// }
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::dumpToParaview(std::string /*filename*/) const {
#ifdef USING_IOHELPER
UInt n = this->size();
const Surface& data = *this;
/* build 3D coordinates array to be usable with iohelper */
// cerr << "allocate coords for paraview" << endl;
Real* coords = new Real[n * n * 3];
// Real * coords = new Real[(2*n-1)*(2*n-1)*3];
// cerr << "allocate zcoords for paraview" << endl;
Real* zcoords = new Real[n * n * 3];
// Real * zcoords = new Real[(2*n-1)*(2*n-1)*3];
// cerr << "allocate zcoords2 for paraview" << endl;
Real* zcoords2 = new Real[n * n * 3];
// Real * zcoords2 = new Real[(2*n-1)*(2*n-1)*3];
int index = 0;
// for (int i=-n+1; i<n; i++) {
// for (int j=-n+1; j<n; j++) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
coords[3 * index] = 1. * i / n;
coords[3 * index + 1] = 1. * j / n;
coords[3 * index + 2] = 0;
zcoords[3 * index] = 0;
zcoords[3 * index + 1] = 0;
zcoords2[3 * index] = 0;
zcoords2[3 * index + 1] = 0;
int i1 = (i % n);
int j1 = (j % n);
if (i1 < 0)
i1 += n;
if (j1 < 0)
j1 += n;
zcoords[3 * index + 2] = data(i1, j1).real();
zcoords2[3 * index + 2] = data(i1, j1).imag();
++index;
}
}
iohelper::DumperParaview dumper;
dumper.setPoints(coords, 3, index, filename.c_str());
dumper.setConnectivity(nullptr, iohelper::POINT_SET, 1, iohelper::FORTRAN_MODE);
dumper.addNodeDataField("data.re", zcoords, 3, index);
dumper.addNodeDataField("data.im", zcoords2, 3, index);
dumper.setPrefix("./paraview");
dumper.init();
dumper.dump();
delete[] coords;
delete[] zcoords;
delete[] zcoords2;
#endif // USING_IOHELPER
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::recenterHeights() {
UInt n = this->size();
Real zmean = SurfaceStatistics::computeAverage(*this);
for (UInt i = 0; i < n * n; ++i)
(*this)(i) = (*this)(i)-zmean;
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Surface<T>::generateWhiteNoiseSurface(Surface<Real>& sources,
long random_seed) {
// set the random gaussian generator
random_engine gen(random_seed);
normal_distribution<Real> gaussian_generator;
UInt n = sources.size();
for (UInt i = 0; i < n; ++i)
for (UInt j = 0; j < n; ++j) {
sources(i, j) = gaussian_generator(gen);
}
}
/* -------------------------------------------------------------------------- */
template class Surface<Real>;
template class Surface<UInt>;
template class Surface<Int>;
template class Surface<bool>;
} // namespace tamaas
c4science · Help