Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F70896288
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
Mon, Jul 8, 03:52
Size
14 KB
Mime Type
text/x-c++
Expires
Wed, Jul 10, 03:52 (2 d)
Engine
blob
Format
Raw Data
Handle
18883801
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