Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/surface_generator_voss.cpp
562a3d65fd33features/netcdf-mpi

Last Change
Raw File

surface_generator_voss.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_voss.hh"
#include "math.h"
#include <iostream>
#include <fstream>
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
#define TWOPOW(A) (pow(2.0,((Real)(A))))
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
using namespace std;
/* -------------------------------------------------------------------------- */
SurfaceGeneratorVoss::SurfaceGeneratorVoss():
SurfaceGenerator(){
}
/* -------------------------------------------------------------------------- */
/* Check if a number is a power of two */
int isPowerOfTwo (unsigned int x)
{
return ((x != 0) && ((x & (~x + 1)) == x));
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorVoss::Init(){
SurfaceGenerator::Init();
int gridsize = surface->size();
/* Post-process Hurst exponent */
if ((Hurst<0) || (Hurst>1)) {
SURFACE_FATAL("Hurst exponent should be between 0.0 and 1.0");
}
if (!isPowerOfTwo(gridsize))
SURFACE_FATAL("Grid size should be a power of two and is " << gridsize);
Nloops = (int)(log(gridsize)/log(2.0));
/* Post-process RMS */
if ((rms < 0) || (fabs(rms) < 1e-10)) {
SURFACE_FATAL("RMS should be larger than zero");
}
/* Initialize specific lists for the surface algorithm */
Nsquares = (int *)malloc((gridsize+1)*sizeof(int));
Ndiamonds = (int *)malloc((gridsize+1)*sizeof(int));
return;
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorVoss::createRoughSurface(Real *err) {
int i, j, k, p, n;
int nsquares, ndiamonds;
int dns, dnd;
int midx, midy;
int top, bot, lft, rgt;
Real rval, roughparam;
Real curval, newval;
Real min, max, avg, scaling;
/* Set some variables */
int gridsize = surface->size();
n = gridsize+1;
roughparam = 1.0;
/* Initiate the dynamic arrays */
for (i=0; i<n; i++) { Nsquares[i]=0; Ndiamonds[i]=0; }
/* Fill the surface arrays with zeros */
for (i=0; i<(n*n); i++) {
surface->at(i) = 0.0;
}
/* Put random values on the corner points */
SURFACE_FATAL("AAAAAAAAAAAAAAAA");
// rval = roughness(roughparam);
surface->at(0) = rval;
surface->at((n*n)-1) = rval;
surface->at(n-1) = rval;
surface->at((n*n)-n) = rval;
/* Start the looping */
for (p=0; p<Nloops; p++) {
/* First the squares */
nsquares = (int) TWOPOW(p);
dns = gridsize/nsquares;
Nsquares[0] = 0; /* Should be zero always, just to be sure */
for (i=1; i<(nsquares+1); i++) { Nsquares[i] = Nsquares[i-1]+dns; }
for (i=0; i<nsquares; i++) {
for (j=0; j<nsquares; j++) {
curval = 0;
curval += surface->at(Nsquares[i] + n*Nsquares[j] );
curval += surface->at(Nsquares[i+1] + n*Nsquares[j] );
curval += surface->at(Nsquares[i] + n*Nsquares[j+1]);
curval += surface->at(Nsquares[i+1] + n*Nsquares[j+1]);
SURFACE_FATAL("AAAAAAAAAAAAAAAA");
// newval = 0.25*curval + roughness(roughparam);
midx = Nsquares[i+1]-dns/2;
midy = Nsquares[j+1]-dns/2;
surface->at(midx + n*midy) += newval;
}
}
/* Then the edges of the squares */
ndiamonds = (int) TWOPOW(p+1);
dnd = gridsize/ndiamonds;
Ndiamonds[0] = 0; /* Should be zero always, just to be sure */
for (i=1; i<(ndiamonds+1); i++) { Ndiamonds[i] = Ndiamonds[i-1]+dnd; }
for (i=0, k=0; i<(ndiamonds+1); i++) {
for (j=0; j<(ndiamonds+1); j++, k++) {
if ((k%2)==1) {
top = Ndiamonds[j]+dnd; top += top<0?gridsize:0; top -= top>gridsize?gridsize:0;
bot = Ndiamonds[j]-dnd; bot += bot<0?gridsize:0; bot -= bot>gridsize?gridsize:0;
lft = Ndiamonds[i]+dnd; lft += lft<0?gridsize:0; lft -= lft>gridsize?gridsize:0;
rgt = Ndiamonds[i]-dnd; rgt += rgt<0?gridsize:0; rgt -= rgt>gridsize?gridsize:0;
curval = 0;
curval += surface->at(rgt + n*(Ndiamonds[j]));
curval += surface->at(lft + n*(Ndiamonds[j]));
curval += surface->at(Ndiamonds[i] + n*(bot) );
curval += surface->at(Ndiamonds[i] + n*(top) );
SURFACE_FATAL("AAAAAAAAAAAAAAAA");
// newval = 0.25*curval + roughness(roughparam);
surface->at(Ndiamonds[i] + n*Ndiamonds[j]) += newval;
}
}
}
/* And now keep things periodic */
for (i=0; i<n; i++) {
surface->at(i+n*(n-1)) = surface->at(i);
surface->at((n-1)+n*i) = surface->at(n*i);
}
/* Update roughness parameter */
roughparam *= TWOPOW(-Hurst);
}
/* Reformat the grid */
min = SurfaceStatistics::computeMinimum(*surface);
max = SurfaceStatistics::computeMaximum(*surface);
avg = SurfaceStatistics::computeAverage(*surface);
scaling = rms/(max-min);
for (i=0; i<(n*n); i++) {
surface->at(i) -= avg;
surface->at(i) *= scaling;
}
/* Error estimate */
*err = fabs(surface->at(0) - SurfaceStatistics::computeAverage(*surface));
return;
}
Surface<Real> & SurfaceGeneratorVoss::buildSurface(){
//int i, j, n;
int n;
Real err;
int gridsize = surface->size();
/* Allocate the surface arrays */
n = gridsize+1;
surface = new Surface<Real>(n,1.);
// heights = (Real *)allocate(n*n*sizeof(Real),"heights");
/* Build the top surface */
fprintf(stdout,"Building surface...\n"); fflush(stdout);
err = rms;
while (err>(rms/4)) { createRoughSurface(&err); }
return *surface;
}
__END_TAMAAS__
c4science · Help