Commit 56d49edb authored by Lukas Riedel's avatar Lukas Riedel

removed old RFG files

parent ce1c4174
/*
* File: FieldData.hh
* Author: A. Ngo (09/2014)
*
* Data structure used to read the input data for the random field to be generated.
* Corresponds to "yfield.cfg".
*
*/
#ifndef FIELD_DATA_HH
#define FIELD_DATA_HH
namespace Dune {
namespace Dorie {
template<typename RF, const int dim>
class FieldData {
public:
Vector<RF> extensions;
Vector<unsigned int> nCells;
Vector<RF> gridsizes;
Vector<RF> correlations;
RF beta;
RF fieldVariance;
std::string variogramModel;
RF embeddingFactor;
bool newField;
bool newEV;
bool showEV;
int seed;
std::string location;
FieldData(){
extensions.resize(dim);
nCells.resize(dim);
gridsizes.resize(dim);
correlations.resize(dim);
};
void read(const Dune::ParameterTree& configuration){
Dune::FieldVector<RF,dim> extensions_ = configuration.get<Dune::FieldVector<RF,dim>>("generator.extensions");
extensions[0] = extensions_[0];
extensions[1] = extensions_[1];
if( dim == 3 )
extensions[2] = extensions_[2];
fieldVariance = configuration.get<RF>("generator.variance");
Dune::FieldVector<RF,dim> nCells_ = configuration.get<Dune::FieldVector<RF,dim>>("generator.fft.N");
nCells[0] = nCells_[0];
nCells[1] = nCells_[1];
if( dim == 3 )
nCells[2] = nCells_[2];
for(int i=0;i<dim;++i)
gridsizes[i] = extensions[i] / RF( nCells[i] );
Dune::FieldVector<RF,dim> correlations_ = configuration.get<Dune::FieldVector<RF,dim>>("generator.fft.correlationLengths");
correlations[0] = correlations_[0];
correlations[1] = correlations_[1];
if( dim == 3 )
correlations[2] = correlations_[2];
beta = configuration.get<RF>("generator.fft.beta");
embeddingFactor = configuration.get<RF>("generator.fft.embeddingFactor");
variogramModel = configuration.get<std::string>("generator.fft.variogramModel");
newField = ("true"==configuration.get<std::string>("generator.fft.newField"));
newEV = ("true"==configuration.get<std::string>("generator.fft.newEV"));
showEV = ("true"==configuration.get<std::string>("generator.fft.showEV"));
seed = configuration.get<int>("generator.fft.seed");
location = configuration.get<std::string>("generator.fft.outputPath");
userInputCheck();
};
void userInputCheck(){
RF nAllCells = nCells[0] * nCells[1];
if(dim == 3)
nAllCells *= nCells[2];
if( nCells[0] <= 0 || nCells[1] <= 0 || nAllCells <= 0 )
DUNE_THROW(IOError,"Random field extensions (N) need to be larger than zero!");
if(nAllCells > 2.5e8)
DUNE_THROW(IOError,"Random fields with more than 2.5e8 cells are not supported!");
//For more information look up https://www.hdfgroup.org/hdf5-quest.html#p2gb
if(correlations[0] <= .0 || correlations[1] <= .0 || (dim == 3 && correlations[2] <= .0))
DUNE_THROW(IOError,"Random field correlation lengths need to be larger than zero!");
};
void printInfos(){
std::cout << "Dimension : " << dim << std::endl;
std::cout << "Domain : " << Vector2CSV( extensions ) << std::endl;
std::cout << "Grid : " << Vector2CSV( nCells ) << std::endl;
std::cout << "cell-size : " << Vector2CSV( gridsizes ) << std::endl;
std::cout << "Variogram model : " << variogramModel << std::endl;
std::cout << "Correlation-lengths : " << Vector2CSV( correlations ) << std::endl;
std::cout << "Field trend : " << beta << std::endl;
std::cout << "Field variance : " << fieldVariance << std::endl;
std::cout << "embedding factor : " << embeddingFactor << std::endl;
std::cout << "seed : " << seed << std::endl;
};
};
}
}
#endif /* FIELD_DATA_HH */
This diff is collapsed.
This diff is collapsed.
/*
* File: Vector.hh
* Author: ngo
*
* Created on July 15, 2010, 2:32 PM
*/
#ifndef _VECTOR_HH
#define _VECTOR_HH
#include <vector>
#include <cmath>
#include <assert.h>
#include <iostream>
#include <sstream>
#include <iomanip>
template<typename ComponentType>
class Vector : public std::vector<ComponentType> {
typedef std::vector<ComponentType> basetype;
public:
Vector() : std::vector<ComponentType>() {
}
Vector( const std::vector<ComponentType>&x )
: std::vector<ComponentType>( x ) {
}
Vector( const Vector& orig ) : basetype(orig) {
}
Vector(const int size, const ComponentType defaultvalue_ = 0)
: std::vector<ComponentType>( size, defaultvalue_ ) {
}
/*
Vector(const int dim, const ComponentType* array_){
for(size_t i=0; i<dim; i++){
ComponentType a = array_[i];
this->push_back(a);
}
}
*/
void reset( const int dim ){
Vector &self = *this;
for( int i=0;i<dim;++i )
self[i] = 0;
}
// inner product with another vector:
ComponentType operator*(Vector & x) {
assert( x.size() == this->size() );
ComponentType sum( 0 );
Vector & self = *this;
for( size_t i = 0; i < this->size(); ++i )
sum += self[i] * x[i];
return sum;
}
// get the product of all the components
ComponentType componentsproduct() const {
ComponentType product( 1 );
const Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
product *= self[i];
return product;
}
// square of the euclidean norm:
ComponentType two_norm_sqr() const {
ComponentType sum( 0 );
const Vector & self = *this;
for (int i = 0; i < (int) this->size(); ++i)
sum += self[i] * self[i];
return sum;
}
// euclidean norm:
ComponentType two_norm() const {
return sqrt(two_norm_sqr());
}
// Manhattan norm:
ComponentType one_norm() const {
ComponentType sum( 0 );
const Vector & self = *this;
for (int i = 0; i < (int) this->size(); ++i)
sum += std::abs( self[i] );
return sum;
}
// Maximum norm:
ComponentType maximum_norm() const {
ComponentType maximum( 0 );
const Vector & self = *this;
for (int i = 0; i < (int) this->size(); ++i)
maximum = std::max( std::abs( self[i] ), maximum );
return maximum;
}
// Arithmetic mean:
ComponentType mean() const {
ComponentType sum( 0 );
const Vector & self = *this;
int N = (int) this->size();
for (int i = 0; i<N; ++i)
sum += self[i];
return REAL(sum) / REAL(N);
}
// Variance:
void mean_and_variance( REAL& mean_value, REAL& variance) const {
ComponentType sum( 0 );
const Vector & self = *this;
mean_value = mean();
int N = (int) this->size();
for (int i = 0; i<N; ++i){
ComponentType summand = self[i] - mean_value;
sum += summand * summand;
}
variance = REAL(sum) / REAL(N-1);
return;
}
// x[i] = value for all components
Vector &operator=(const ComponentType value) {
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] = value;
return *this;
}
// x[i] *= value
Vector &operator*=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] *= value;
return *this;
}
// x[i] /= value
Vector &operator/=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] /= value;
return *this;
}
// x[i] += value
Vector &operator+=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] += value;
return *this;
}
// x += y
Vector &operator+=(const Vector& b) {
Vector& self = *this;
assert( self.size() == b.size() );
for (size_t i = 0; i < this->size(); ++i)
self[i] += b[i];
return *this;
}
// x -= y
Vector &operator-=(const Vector& b) {
Vector& self = *this;
assert( self.size() == b.size() );
for (size_t i = 0; i < this->size(); ++i)
self[i] -= b[i];
return *this;
}
// x[i] -= value
Vector &operator-=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] -= value;
return *this;
}
Vector operator-(const Vector& b){
const Vector &a = *this;
assert( a.size() == b.size() );
Vector c(a);
assert( c.size() == b.size() );
c.axpy(-1,b);
return c;
};
// x = x + alpha * y :
Vector & axpy(const ComponentType alpha, const Vector & y) {
assert( this->size() == y.size());
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] += alpha * y[i];
return *this;
}
// x = x + alpha * y
Vector & scaled_add(const ComponentType alpha, const Vector & y) {
assert( this->size() == y.size());
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] += alpha * y[i];
return *this;
}
// x = alpha * x + y
Vector & a_times_x_plus_y(const ComponentType alpha, const Vector & y) {
assert( this->size() == y.size());
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] = alpha * self[i] + y[i];
return *this;
}
void print(void) const {
const Vector &self = *this;
std::cout << "( ";
for (size_t i = 0; i < self.size(); i++) {
std::cout << self[i] << " ";
}
std::cout << " )" << std::endl;
}
};
template<typename ComponentType>
std::ostream & operator <<(std::ostream & os, const Vector<ComponentType> & x) {
os << "( ";
for (int r = 0; r < (int) x.size(); ++r) {
os << x[r] << " ";
}
os << ")" << std::endl;
return os;
}
#endif /* _VECTOR_HH */
typedef double CTYPE;
typedef double REAL;
typedef unsigned int UINT;
typedef int INT;
template<int dim>
struct P0Layout
{
bool contains (Dune::GeometryType gt)
{
if (gt.dim()==dim) return true;
return false;
}
};
struct FEMType {
enum Type{ DG, CG };
};
// define the data type of the HDF5 files here. 32 or 64 bit!!!!
#define HDF5_DATA_TYPE H5T_IEEE_F64LE //define for 64 bit machine
//#define HDF5_DATA_TYPE H5T_IEEE_F32LE //define for 32 bit machine
This diff is collapsed.
#ifndef DUNE_FUNCEP_FIELDGENTOOLS_HH
#define DUNE_FUNCEP_FIELDGENTOOLS_HH
namespace Dune {
namespace Dorie {
UINT indexconversion_to_1d( const Vector<UINT>& index,
const Vector<UINT>& nCells )
{
// "http://www.fftw.org/doc/Dynamic-Arrays-in-C.html#Dynamic-Arrays-in-C"!
if( nCells.size() > 2 )
return index[0] + nCells[0] * ( index[1] + nCells[1] * index[2] ); // 3d to 1d
else
return index[0] + nCells[0] * index[1]; // 2d to 1d
};
UINT indexconversion_3d_to_1d( const UINT iz
, const UINT iy
, const UINT ix
, const UINT Nz
, const UINT Ny
, const UINT Nx
) {
// "http://www.fftw.org/doc/Dynamic-Arrays-in-C.html#Dynamic-Arrays-in-C"
UINT l = ix + Nx * ( iy + Ny*iz );
return l;
};
// index conversion from 2d-index to 1d-index
UINT indexconversion_2d_to_1d(
const UINT iy
, const UINT ix
, const UINT Ny
, const UINT Nx
)
{
// UINT l = iy * Nx + ix;
UINT l = ix + Nx * iy;
return l;
};
template<typename BVEC, typename VEC>
void copy_to_std( const BVEC& backend_vector, VEC& standard_vector ){
for( auto it = backend_vector.begin(); it!=backend_vector.end(); ++it )
standard_vector.push_back(*it);
}
template<typename BVEC, typename VEC>
void copy_from_std( const VEC& standard_vector, BVEC& backend_vector ){
for( int i=0; i<standard_vector.size(); ++i )
backend_vector.base()[i] = standard_vector[i];
}
template<typename V>
std::string Vector2CSV( const V& v ){
std::stringstream output;
unsigned int j=0;
do{
output << v[j] << ", ";
++j;
} while(j<v.size()-1);
output << v[j];
return output.str();
}
void deleteFile( const std::string filename ){
if( 0 != remove( filename.c_str() ) ){
std::cout << "Error deleting file: " << filename << std::endl;
}
else{
std::cout << "File deleted: " << filename << std::endl;
}
};
/*
* This simple function checks if a file exists.
*/
bool bFileExists( const std::string& filename )
{
struct stat buf;
if( stat( filename.c_str(), &buf ) != -1 ){
UINT filesize = (UINT) buf.st_size;
std::cout << "File found: " << filename
<< ", file size = " << filesize
<< std::endl;
if( filesize > 0 ){
return true;
}
else{
deleteFile( filename );
return false;
}
}
else{
return false;
}
}
} // end namespace FunCEP
} // end namespace Dune
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment