"""
Classes for generating pseudo-random surfaces based on description of FFT:
===========================================================================
===========================================================================
Each class inherits functionality from the Surface but changes the
__init__ and discretise functions
===========================================================================
===========================================================================
DiscFreqSurface:
Generate a surface containing only specific frequency components
ProbFreqSurface:
Generate a surface containing normally distributed amplitudes with a
specified function for the variance of the distribution based on the
frequency
DiscFreqSurface:
Generate a surface containing frequency components with amplitude
specified by a function of the frequency
===========================================================================
===========================================================================
"""
import typing
from numbers import Number
import numpy as np
from .Surface_class import _AnalyticalSurface, Surface
__all__ = ["DiscFreqSurface", "ProbFreqSurface", "HurstFractalSurface"]
[docs]class DiscFreqSurface(_AnalyticalSurface):
r""" Surfaces with discrete frequency components
This surface produces a profile with frequency components which are constant in th y direction,
these can be rotated and combined by addition or subtraction to give any combination of frequencies.
Parameters
----------
frequencies: Sequence[float]
The frequencies present in the surface
amplitudes: Sequence[float], optional ((1, ))
The amplitude of each frequency, must be same length as frequencies
phases: typing.Sequence[float] = (0,)
The phases of each of the frequencies, must be the same length as frequencies
rotation: Number = 0
If set the surface is rotated by the set amount, in radians
shift: tuple, optional (None)
If set the origin of the surface is shifted by the set amount in the x and y directions, should be a two element
tuple of float. By default, the surface is shifted so that the origin becomes the centre of the surface. To stop
this behaviour specify a shift of (0, 0)
generate: bool, optional (False)
If True the surface profile is generated on instantiation, else it can be generated by the discretise method.
If True two of the: grid_spacing, extent or shape must also be set
grid_spacing: float, optional (None)
The grid spacing of the surface profile
extent: tuple, optional (None)
The overall dimensions of the surface in the same units as the grid spacing, should be a two element tuple of
float
shape: tuple = None
The number of points in each direction of the surface array, should be a two element tuple of integers
See Also
--------
ProbFreqSurface
HurstFractalSurface
Notes
-----
Roughness functions are aliased from the functions provided in the surface
module
Examples
--------
>>> mySurf=DiscFreqSurface(10, 0.1)
>>> mySurf.extent=[0.5,0.5]
>>> mySurf.discretise(0.001)
"""
is_discrete = False
surface_type = 'discreteFreq'
[docs] def __init__(self, frequencies: typing.Sequence[float], amplitudes: typing.Sequence[float] = (1,),
phases: typing.Sequence[float] = (0,), rotation: Number = 0,
shift: typing.Optional[tuple] = None,
generate: bool = False, grid_spacing: float = None,
extent: tuple = None, shape: tuple = None):
if type(frequencies) is list or type(frequencies) is np.ndarray:
self.frequencies = frequencies
else:
raise ValueError('Frequencies, amplitudes and phases must be equal'
'length lists or np.arrays')
is_complex = [type(amp) is complex for amp in amplitudes]
if any(is_complex):
if not len(frequencies) == len(amplitudes):
raise ValueError('Frequencies, amplitudes and phases must be'
' equal length lists or np.arrays')
else:
self.amplitudes = amplitudes
else:
if not len(frequencies) == len(amplitudes) == len(phases):
raise ValueError('Frequencies, amplitudes and phases must be'
' equal length lists or np.arrays')
else:
cplx_amps = []
for idx in range(len(amplitudes)):
cplx_amps.append(amplitudes[idx] *
np.exp(1j * phases[idx]))
self.amplitudes = cplx_amps
super().__init__(generate=generate, rotation=rotation, shift=shift,
grid_spacing=grid_spacing, extent=extent, shape=shape)
def _height(self, x_mesh, y_mesh):
profile = np.zeros_like(x_mesh)
for idx in range(len(self.frequencies)):
profile += np.real(self.amplitudes[idx] *
np.exp(-1j * self.frequencies[idx] * x_mesh * 2 * np.pi))
return profile
def __repr__(self):
string = self._repr_helper()
return f'DiscFreqSurface({self.frequencies}, {self.amplitudes}{string})'
[docs]class ProbFreqSurface(_AnalyticalSurface):
"""
ProbFreqSurface(H, qr, qs)
Generates a surface with all possible frequencies in the fft represented
with amplitudes described by the probability distribution given as input.
Defaults to the parameters used in the contact mechanics challenge
This class only works for square 2D domains
For more information the definitions of the input parameters refer to
XXXXXX contact mechanics challenge paper
"""
is_discrete = False
surface_type = 'Random'
[docs] def __init__(self, h=2, qr=0.05, qs=10,
generate: bool = False, grid_spacing: float = None,
extent: tuple = None, shape: tuple = None):
self.h = h
self.qs = qs
self.qr = qr
super().__init__(grid_spacing=grid_spacing, extent=extent, shape=shape, generate=generate)
def rotate(self, radians: Number):
raise NotImplementedError("Probabilistic frequency surface cannot be rotated")
def shift(self, shift: tuple = None):
if shift is None:
return
raise NotImplementedError("Probabilistic frequency surface cannot be shifted")
def _height(self, x_mesh, y_mesh):
grid_spacing, extent, shape = check_coords_are_simple(x_mesh, y_mesh)
qny = np.pi / grid_spacing
u = np.linspace(0, qny, shape[0])
u_mesh, v_mesh = np.meshgrid(u, u)
freqs = np.abs(u_mesh + v_mesh)
varience = np.zeros(freqs.shape)
varience[np.logical_and((1 / freqs) > (1 / self.qr), (2 * np.pi / freqs) <= (extent[0]))] = 1
varience[np.logical_and((1 / freqs) >= (1 / self.qs), (1 / freqs) < (1 / self.qr))] = \
(freqs[np.logical_and(1 / freqs >= 1 / self.qs, 1 / freqs < 1 / self.qr)] / self.qr) ** (-2 * (1 + self.h))
fou_trans = np.reshape(np.array([np.random.normal() * var ** 0.5 for var in varience.flatten()]), freqs.shape)
return np.real(np.fft.ifft2(fou_trans))
def __repr__(self):
string = self._repr_helper()
return f'ProbFreqSurface(h={self.h}, qr={self.qr}, qs={self.qs}{string})'
[docs]class HurstFractalSurface(Surface):
r"""Hurst fractal surfaces
Parameters
----------
sigma: float
The RSM roughness of the surface
hurst_exponent: float
The hurst exponent, must be between 0 and 1, related to the fractal dimension by D = 3-H
roll_off_frequency: float, optional (0.0)
generate: bool, optional (False)
If true the surface profile is generated on instantiation, two of: grid_spacing, extent or shape must be set
grid_spacing: float, optional (None)
The grid spacing of the surface profile
extent: tuple, optional (None)
The overall surface dimensions in the x and y directions
shape: tuple, optional (None)
The number of grid points in the x and y directions, computation is faster for powers of 2
See Also
--------
ProbFreqSurface
RandomSurface
surface_like
Notes
-----
generates a hurst fractal surface with frequency components from q0 to
cut off frequency in even steps of q0.
amplitudes are given by:
q0 amplitude\*\*2 \*((h\*\*2+k\*\*2)/2)\*\*(1-Hurst parameter)
where h,k = \-N...N
where N=cut off frequency/ q0
phases are randomly generated on construction of the surface object,
repeated calls to the discretise function will discretise on the same surface
but repeated calls to this class will generate different realisations
References
----------
A new efficient numerical method for contact mechanics of rough surfaces
C.Putignano L.Afferrante G.Carbone G.Demelio
Examples
--------
>>> #create the surface object with the specified fractal parameters
>>> my_surface=HurstFractalSurface(1,0.2,1000, shape=(128, 128), grid_spacing=0.01)
>>> #descrtise the surface over a grid 1 unit by 1 unit with a grid_spacing of 0.01
>>> my_surface.discretise()
>>> my_surface.show()
"""
is_discrete = False
surface_type = "hurstFractal"
[docs] def __init__(self, sigma: float, hurst_exponent: float, roll_off_frequency: float = 0, generate: bool = False,
grid_spacing: float = None, extent: tuple = None, shape: tuple = None):
self.input_params = (sigma, hurst_exponent, roll_off_frequency)
if hurst_exponent > 1 or hurst_exponent < 0:
raise ValueError('Hurst exponent must be between 0 and 1')
self._hurst_exponent = hurst_exponent
self._sigma = sigma
self._roll_off_frequency = roll_off_frequency
super().__init__(grid_spacing=grid_spacing, extent=extent, shape=shape)
if generate:
self.discretise()
def discretise(self, return_new: bool = False):
"""Generate a new profile realisation, return a new surface if needed
Parameters
----------
return_new: bool
If True a new surface instance is returned, else the profile property of the current surface is set
Returns
-------
out: Surface
A new surface object with the profile set, will have the same shape and grid spacing as the current object
only returned if return_new is True
Notes
-----
As a side effect this will set the FFT and PSD properties of the discretised surface.
Copyright (c) 2016, Mona Mahboob Kanafi
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution
* Neither the name of Aalto University nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
References
----------
Code ported from the matlab version:
https://uk.mathworks.com/matlabcentral/fileexchange/60817-surface-generator-artificial-randomly-rough-surfaces?s_tid=mwa_osa_a
"""
if self.profile is not None and not return_new:
raise ValueError('Profile is already set, set the return_new argument to true to return a new surface '
'instance with the same fractal properties')
if self.shape is None or self.grid_spacing is None:
raise ValueError('Grid spacing and shape of the surface must be set before a hurst fractal can be '
'generated')
m, n = [s + s % 2 for s in self.shape]
grid_spacing = self.grid_spacing
hurst_exponent = self._hurst_exponent
sigma = self._sigma
qr = self._roll_off_frequency
lx, ly = [n_pts * grid_spacing for n_pts in [m, n]]
# make the wave vectors
qx = np.array([2 * np.pi / m * k for k in range(m)])
qx = np.unwrap((np.fft.fftshift(qx)) - 2 * np.pi) / grid_spacing
qy = np.array([2 * np.pi / n * k for k in range(n)])
qy = np.unwrap((np.fft.fftshift(qy)) - 2 * np.pi) / grid_spacing
# find absolute frequencies
qx_mesh, qy_mesh = np.meshgrid(qx, qy)
rho = (qy_mesh ** 2 + qx_mesh ** 2) ** 0.5
# make the PSD matrix
psd = rho ** (-2 * (hurst_exponent + 1))
psd[rho < qr] = qr ** (-2 * (hurst_exponent + 1)) if qr > 0 else np.inf
psd[n // 2, m // 2] = 0.0
# Scale surface to sigma (apply rms)
rms_f2d = np.sqrt(sum(psd.flatten()) * (2 * np.pi) ** 2 / lx / ly)
alpha = sigma / rms_f2d
psd *= alpha ** 2
self.psd = np.fft.ifftshift(psd)
# converting the PSD to the fft
magnitudes = np.sqrt(psd / (grid_spacing ** 2 / (n * m * (2 * np.pi) ** 2)))
magnitudes = _conj_sym(magnitudes)
phases = -np.pi + (2 * np.pi) * np.random.rand(n, m)
phases = _conj_sym(phases, neg=True)
u, v = _pol2cart(magnitudes, phases)
fou_trans = np.empty(u.shape, dtype=complex)
fou_trans.real = u
fou_trans.imag = v
fou_trans = np.fft.ifftshift(fou_trans)
profile = np.fft.ifft2(fou_trans)
profile = profile.real
if return_new:
out = Surface(profile=profile, grid_spacing=grid_spacing)
out.psd = psd
out.fft = fou_trans
return out
self.profile = profile
self.fft = fou_trans
def __repr__(self):
string = super().__repr__()[:-1]
input_param_string = ', '.join([str(p) for p in self.input_params])
return f'HurstFractal{string}, {input_param_string[:-1]})'
def rotate(self, radians: Number):
raise NotImplementedError("Hurst fractal surface cannot be rotated")
def shift(self, shift: tuple = None):
if shift is None:
return
raise NotImplementedError("Hurst fractal surface cannot be shifted")
def check_coords_are_simple(x_mesh, y_mesh):
"""
Checks that the coordinates are of the type np.meshgrid(np.arrange(0, end, step), np.arrange(0, end, step))
Parameters
----------
x_mesh: np.ndarray
y_mesh: np.ndarray
Returns
-------
grid_spacing, extent, shape
"""
x_mesh_check, y_mesh_check = np.meshgrid(x_mesh[0, :], x_mesh[0, :])
difference = np.diff(x_mesh[0, :])
if not np.allclose(difference, difference[0]):
raise ValueError('x and y points must be evenly spaced and sorted')
if not x_mesh[0, 0] == 0:
raise ValueError('x and y points must start at 0')
if not x_mesh_check == x_mesh and y_mesh_check == y_mesh:
raise ValueError("For discretising a Probabilistic frequency surface the x and y coordinates should form"
" an evenly spaced grid aligned with the axes")
extent = (x_mesh[0, -1], x_mesh[0, -1])
grid_spacing = difference[0]
shape = [ex / grid_spacing + 1 for ex in extent]
return grid_spacing, extent, shape
def _conj_sym(matrix: np.ndarray, in_place=True, neg=False):
""" apply conjugate symmetry to a matrix
"""
try:
n, m = matrix.shape
except AttributeError:
raise ValueError(f'Input matrix must be a 2d numpy array, got {type(matrix)}')
mult = -1 if neg else 1
if not in_place:
matrix = matrix.copy()
matrix[0, 0] = 0.0
matrix[0, m // 2] = 0.0
matrix[n // 2, 0] = 0.0
matrix[n // 2, m // 2] = 0.0
matrix[1:, 1:m // 2] = mult * np.rot90(matrix[1:, m // 2 + 1:], 2)
matrix[0, 1:m // 2] = mult * np.flip(matrix[0, m // 2 + 1:]) # just needs to flip it
matrix[n // 2 + 1:, 0] = mult * np.flip(matrix[1:n // 2, 0]) # also just needs to flip it
matrix[n // 2 + 1:, m // 2] = mult * np.flip(matrix[1:n // 2, m // 2])
return matrix
def _pol2cart(rho, phi):
x = rho * np.cos(phi)
y = rho * np.sin(phi)
return x, y