Homec4science
Diffusion LibMultiScale (features/cadd) Browse

libmultiscale/src/geometry/surface.hh
88abcabc76bafeatures/cadd

Last Change
Raw File

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__ */
c4science · Help