Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/surface_generator_filter_fft.cpp
7d473b01cbbafeatures/declarative

Last Change
Raw File

surface_generator_filter_fft.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 "surface_generator_filter_fft.hh"
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <cmath>
#include <fstream>
#include <fftw3.h>
/* -------------------------------------------------------------------------- */
SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT():
SurfaceGeneratorFilter(){
}
/* -------------------------------------------------------------------------- */
complex SurfaceGeneratorFilterFFT::computeFractalPSDProfile(int i, int j){
Real q = sqrt(i*i+j*j);
complex value(Real(0.),Real(0.));
if (q < q0) value = 0.;
else if (q <= q1) value = 1.;
else if (q > q2) value = 0.;
else value = pow(q/q1,-2.*(hurst+1));
return value;
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorFilterFFT::buildCorrelationInput(){
int n = surface->size();
// all entries are set to zero automatically
Surface<Real> power_input(n,1.);
std::cout << "rms: " << rms << " , hurst: " << hurst << " , q0: " << q0 << " , q1: " << q1 << " , q2: " << q2 << std::endl;
power_input(0,0) = 0.;
for (int i = 1 ; i < n/2 ; ++i)
for (int j = 1 ; j < n/2 ; ++j){
complex value = computeFractalPSDProfile(i,j);
power_input(i,j) = value.real();
power_input(n-i,n-j) = value.real();
power_input(n-i,j) = value.real();
power_input(i,n-j) = value.real();
}
// external line (i or j == 0)
for (int i = 1 ; i < n/2 ; ++i){
complex value = computeFractalPSDProfile(i,0);
power_input(i,0) = value.real();
power_input(0,i) = value.real();
power_input(0,n-i) = value.real();
power_input(n-i,0) = value.real();
value = computeFractalPSDProfile(i,n/2);
power_input(i,n/2) = value.real();
power_input(n-i,n/2) = value.real();
power_input(n/2,i) = value.real();
power_input(n/2,n-i) = value.real();
}
{
power_input(n/2,n/2) = computeFractalPSDProfile(n/2,n/2).real();
power_input(0,n/2) = computeFractalPSDProfile(0,n/2).real();
power_input(n/2,0) = computeFractalPSDProfile(n/2,0).real();
}
power_input.dumpToParaview("power_spectrum_input");
Real average_spectrum = SurfaceStatistics::computeAverageAmplitude(*white_noise_FFT);
std::cout << "white noise average spectrum is " << average_spectrum << std::endl;
// compute fourier transform of filter coefficients
h_coeff_power = new SurfaceComplex<Real>(n,1.);
Real C = average_spectrum;
for (int i = 0 ; i < n ; ++i)
for (int j = 0 ; j < n ; ++j){
h_coeff_power->at(i,j) = sqrt(power_input(i,j)/C);
}
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorFilterFFT::computeFilterCoefficients(){
buildCorrelationInput();
}
/* -------------------------------------------------------------------------- */
c4science · Help