Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80388523
surface.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 31, 11:38
Size
14 KB
Mime Type
text/x-c++
Expires
Mon, Sep 2, 11:38 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20358787
Attached To
rLIBMULTISCALE LibMultiScale
surface.hh
View Options
/**
* @file surface.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Thu Feb 06 09:26:39 2014
*
* @brief Pointwise defined surface
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_MYSURFACE_HH__
#define __LIBMULTISCALE_MYSURFACE_HH__
/* -------------------------------------------------------------------------- */
#include "math_tools.hh"
#include <fstream>
#include <iostream>
#include <math.h>
#include <sstream>
#include <string>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class MySurface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MySurface();
MySurface(UInt _size);
~MySurface();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
//! load heights array from text file usable with plot programs
void loadSurfaceFromTextFile(const std::string &filename,
bool extend_for_periodicity = false);
//! rescale heights for having peak2peak requested
void enforcePeak2Peak(Real p2p_req);
//! rescale heights to enforce height RMS
void enforceHeightRMS(Real rms_req);
//! rescale heights to enforce slope RMS
void enforceSlopeRMS(Real rms_req);
//! rescale heights
void rescaleHeights(Real factor);
//! plot surface in a text file
void plotSurface();
//! detect min and max heights
void computeHeightMinMax();
//! compute average heights of the surface
void computeAverageHeight();
//! compute the RMS of slopes for the loaded surface
void computeRMSSlope();
//! compute the RMS of heights for the loaded surface
void computeRMSHeight();
//! return true if corrdinates provided are under the plane surface (in Z
//! direction)
template <UInt Dim> bool isUnderSurface(const Vector<Dim> &X);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
//! return maximal maplitude of the height profile
Real getAmplitudeMax();
//! return average of heights
Real getAverageHeight();
//! return RMS of slopes
Real getRMSSlope();
//! return RMS of heights
Real getRMSHeight();
//! set scale parameter
template <UInt Dim> void setScale(const Vector<Dim> &_scale);
//! set scale parameter by passing length of the domain
void setLength(Real _length);
//! set grid size and allocate the height array
void setGridSize(UInt _size);
//! change the heights array
void setHeights(std::vector<Real> &_heights);
private:
bool isUnderPlane(Real *p1, Real *p2, Real *p3, Real *p4, Real *p);
bool isUnderLine(Real *p1, Real *p2, Real *p);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
//! scaling factor from grid point
Real scale[2];
//! maximal height value
Real zmax;
//! minimal height value
Real zmin;
//! array of heights
std::vector<Real> heights;
//! grid size
UInt size;
//! rms of slopes
Real rms_slopes;
//! average of heights
Real average_heights;
//! rms of heights
Real rms_heights;
};
/* -------------------------------------------------------------------------- */
inline void MySurface::loadSurfaceFromTextFile(const std::string &filename,
bool extend_for_periodicity) {
std::string line;
std::ifstream myfile(filename.c_str());
size = 0;
if (myfile.is_open()) {
while (!myfile.eof()) {
getline(myfile, line);
std::stringstream s;
s << line;
UInt i, j;
s >> i;
s >> j;
if (size < i)
size = i;
if (size < j)
size = j;
}
++size;
DUMPBYPROC("read " << size << " x " << size << " grid points", DBG_MESSAGE,
0);
if (extend_for_periodicity) {
DUMPBYPROC("extend the grid to create a periodic boundary surface",
DBG_MESSAGE, 0);
++size;
}
myfile.clear(); // forget we hit the end of file
myfile.seekg(0, std::ios::beg); // move to the start of the file
setGridSize(size);
while (!myfile.eof()) {
getline(myfile, line);
std::stringstream s;
s << line;
UInt i, j;
Real h;
s >> i;
s >> j;
s >> h;
heights[i + j * size] = h;
}
myfile.close();
// if need to extend reloop over the edges points
if (extend_for_periodicity) {
// the corners
heights[size - 1] = heights[size * (size - 1)] =
heights[size * size - 1] = heights[0];
// first edge to be copied
for (UInt i = 1; i < size - 1; ++i) {
heights[size * (size - 1) + i] = heights[i];
}
// second edge to be copied
for (UInt i = 1; i < size - 1; ++i) {
heights[i * size + size - 1] = heights[i * size];
}
}
} else
LM_FATAL("cannot open file " << filename);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline bool MySurface::isUnderSurface(const Vector<Dim> &X) {
Real Xp[3];
for (UInt i = 0; i < 2; ++i)
Xp[i] = X[i] / scale[i] * (size - 1);
Xp[2] = X[2];
DUMP("here are the Xps " << Xp[0] << " " << Xp[1] << " " << Xp[2],
DBG_DETAIL);
UInt Xp_upper[2];
UInt Xp_lower[2];
for (UInt i = 0; i < 2; ++i) {
Xp_upper[i] = (int)ceil(Xp[i]);
Xp_lower[i] = (int)floor(Xp[i]);
if (Xp_upper[i] >= size || Xp_lower[i] >= size)
LM_FATAL("this should not happen "
<< X[0] << " " << X[1] << " " << Xp_lower[0] << " "
<< Xp_lower[1] << " " << Xp_upper[0] << " " << Xp_upper[1]);
}
Real local_zmin = std::min(heights[Xp_lower[0] + size * Xp_upper[1]],
heights[Xp_upper[0] + size * Xp_upper[1]]);
local_zmin = std::min(local_zmin, heights[Xp_upper[0] + size * Xp_lower[1]]);
local_zmin = std::min(local_zmin, heights[Xp_lower[0] + size * Xp_lower[1]]);
Real local_zmax = std::max(heights[Xp_lower[0] + size * Xp_upper[1]],
heights[Xp_upper[0] + size * Xp_upper[1]]);
local_zmax = std::max(local_zmax, heights[Xp_upper[0] + size * Xp_lower[1]]);
local_zmax = std::max(local_zmax, heights[Xp_lower[0] + size * Xp_lower[1]]);
if (Xp[2] <= local_zmin)
return true;
if (Xp[2] > local_zmax)
return false;
if (Xp_upper[0] == Xp_lower[0]) {
Real p[2] = {Xp[1], Xp[2]};
Real p1[2] = {Real(Xp_upper[1]), heights[Xp_upper[0] + size * Xp_upper[1]]};
Real p2[2] = {Real(Xp_lower[1]), heights[Xp_upper[0] + size * Xp_lower[1]]};
return isUnderLine(p1, p2, p);
} else if (Xp_lower[1] == Xp_upper[1]) {
Real p[2] = {Xp[0], Xp[2]};
Real p1[2] = {Real(Xp_upper[0]), heights[Xp_upper[0] + size * Xp_upper[1]]};
Real p2[2] = {Real(Xp_lower[0]), heights[Xp_lower[0] + size * Xp_upper[1]]};
return isUnderLine(p1, p2, p);
}
bool test;
Real p[3] = {Xp[0], Xp[1], Xp[2]};
{
Real p1[3] = {Real(Xp_lower[0]), Real(Xp_lower[1]),
heights[Xp_lower[0] + size * Xp_lower[1]]};
Real p2[3] = {Real(Xp_upper[0]), Real(Xp_lower[1]),
heights[Xp_upper[0] + size * Xp_lower[1]]};
Real p3[3] = {Real(Xp_lower[0]), Real(Xp_upper[1]),
heights[Xp_lower[0] + size * Xp_upper[1]]};
Real p4[3] = {Real(Xp_upper[0]), Real(Xp_upper[1]),
heights[Xp_upper[0] + size * Xp_upper[1]]};
test = isUnderPlane(p1, p2, p3, p4, p);
}
return test;
}
/* -------------------------------------------------------------------------- */
inline bool MySurface::isUnderLine(Real *p1, Real *p2, Real *p) {
// write p1, p2 in form of a*p[0] +b*p[1] +c = 0;
Real a = p2[1] - p1[1];
Real b = p1[0] - p2[0];
Real c = p2[0] * p1[1] - p1[0] * p2[1];
if (a * p[0] + b * p[1] + c < 0)
return false;
return true;
}
/* -------------------------------------------------------------------------- */
inline bool MySurface::isUnderPlane(Real *p1, Real *p2, Real *p3, Real *p4,
Real *p) {
// write p1, p2, p3 in form of a*p[0] + b*p[1] + c*p[2] + d = 0;
Real x = p[0] - p1[0];
Real y = p[1] - p2[1];
Real z = p[2];
Real N1 = (-1 + y) * (-1 + x);
Real N2 = -x * (-1 + y);
Real N3 = -y * (-1 + x);
Real N4 = x * y;
Real z_interpolated = N1 * p1[2] + N2 * p2[2] + N3 * p3[2] + N4 * p4[2];
if (z_interpolated < z)
return false;
return true;
}
/* -------------------------------------------------------------------------- */
inline void MySurface::computeHeightMinMax() {
zmax = -1.0;
zmin = 1.0;
for (UInt i = 0; i < size; i++)
for (UInt j = 0; j < size; j++) {
zmax = std::max(heights[i + size * j], zmax);
zmin = std::min(heights[i + size * j], zmin);
}
}
/* -------------------------------------------------------------------------- */
inline void MySurface::computeRMSSlope() {
UInt i, j, k, m, n;
Real *x = &heights[0];
UInt xl, xh, yl, yh;
Real avg;
Real gx, gy;
UInt nn = size * size;
n = size;
// if ((n*n)!=nn) my_error("Something wrong here...");
std::vector<Real> grad(nn);
for (i = 0, m = 0; i < n; i++) {
for (j = 0; j < n; j++) {
k = i + n * j;
xl = ((i == 0) ? (k + n - 1) : (k - 1));
xh = ((i == (n - 1)) ? (k - n + 1) : (k + 1));
yl = ((j == 0) ? (k + nn - n) : (k - n));
yh = ((j == (n - 1)) ? (k - nn + n) : (k + n));
// loc = (*(x+k));
// tmp = ((*(x+xl))-loc);
// tmp += ((*(x+xh))-loc);
// tmp += ((*(x+yl))-loc);
// tmp += ((*(x+yh))-loc);
// slp += (tmp*tmp)/4;
gx = (*(x + xh) - *(x + xl)) / 2;
gy = (*(x + yh) - *(x + yl)) / 2;
grad[m++] = (gx * gx) + (gy * gy);
}
}
avg = MathTools::average(grad);
rms_slopes = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */
// slp = stdev(grad,avg,nn);
}
/* -------------------------------------------------------------------------- */
inline void MySurface::computeAverageHeight() {
average_heights = MathTools::average(heights);
}
/* -------------------------------------------------------------------------- */
inline void MySurface::computeRMSHeight() {
rms_heights = MathTools::stdev(heights, average_heights);
}
/* -------------------------------------------------------------------------- */
inline MySurface::MySurface() {
for (UInt i = 0; i < 2; ++i)
scale[i] = NAN;
zmax = NAN;
;
zmin = NAN;
;
size = UINT_MAX;
rms_slopes = NAN;
;
average_heights = NAN;
;
rms_heights = NAN;
;
}
/* -------------------------------------------------------------------------- */
inline MySurface::MySurface(UInt _size) { setGridSize(_size); }
/* -------------------------------------------------------------------------- */
inline MySurface::~MySurface() {}
/* -------------------------------------------------------------------------- */
inline void MySurface::enforcePeak2Peak(Real p2p_req) {
computeHeightMinMax();
Real actual_p2p = zmax - zmin;
if (actual_p2p == 0)
return;
Real factor = p2p_req / actual_p2p;
rescaleHeights(factor);
computeHeightMinMax();
}
/* -------------------------------------------------------------------------- */
inline void MySurface::enforceHeightRMS(Real rms_req) {
computeAverageHeight();
computeRMSHeight();
if (rms_heights == 0)
return;
Real factor = rms_req / rms_heights;
rescaleHeights(factor);
}
/* -------------------------------------------------------------------------- */
inline void MySurface::enforceSlopeRMS(Real rms_req) {
computeRMSSlope();
if (rms_slopes == 0)
return;
Real factor = rms_req / rms_slopes;
rescaleHeights(factor);
}
/* -------------------------------------------------------------------------- */
inline void MySurface::rescaleHeights(Real factor) {
for (UInt i = 0; i < size * size; ++i) {
heights[i] *= factor;
}
}
/* -------------------------------------------------------------------------- */
inline Real MySurface::getAmplitudeMax() { return (zmax - zmin); }
/* -------------------------------------------------------------------------- */
inline Real MySurface::getAverageHeight() { return average_heights; }
/* -------------------------------------------------------------------------- */
inline Real MySurface::getRMSSlope() { return rms_slopes; }
/* -------------------------------------------------------------------------- */
inline Real MySurface::getRMSHeight() { return rms_heights; }
/* -------------------------------------------------------------------------- */
template <UInt Dim>
inline void MySurface::setScale(const Vector<Dim> &_scale) {
for (UInt i = 0; i < 2; ++i) {
scale[i] = _scale[i];
}
}
/* -------------------------------------------------------------------------- */
inline void MySurface::setGridSize(UInt _size) {
size = _size;
heights.resize(size * size);
}
/* -------------------------------------------------------------------------- */
inline void MySurface::setHeights(std::vector<Real> &_heights) {
heights = _heights;
}
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_MYSURFACE_HH__ */
Event Timeline
Log In to Comment