ContactMechanics package

Subpackages

Submodules

ContactMechanics.FFTElasticHalfSpace module

Implement the FFT-based elasticity solver of ContactMechanics

Convention used for the DFT :

In addition to the sum of the product with the exponential function, the one has to divide by \(n_x n_y\) once during the roundtrip.

When this is actually done is arbitrary.

Our convension:

fourier transform:

\[\tilde h_{op} = \sum_{mn} h_{mn} e^{-i x_{mn} q_{op}}\]

corresponding np.fft.rfft and fftengine.fft

fourier space input fields are assumed to be linked to the realspace field through this fourier transform.

fourier inverse transform:

\[h_{mn} = \frac{1}{n_x n_y} \sum_{op} \tilde h_{op} e^{i x_{mn} q_{op}}\]

corresponding np.fft.irfft and fftengine.fft * fftengine.normalisation

Note that this is different from the definition in Jacobs, T. D. B. et al. Surf. Topogr.: Metrol. Prop. 5, 013001 (2017) (Equations A.3, A.4), that is closer to the continuous fourier transform.

Parseval’s theorem, Convolutions and powers:

The prefactors in front of the sums depend on the definition of the fourier transform.

Convolution theorem:

\[(x * y)_m = \sum_n x_n y_{m-n} = IDFT(\tilde x_k \tilde y_k)_m\]

The power theorem can be deduced from the convolution theorem and states that:

\[\sum_n x_n \overline{y_n} = \frac{1}{N} \sum_n \tilde x_n \overline{\tilde y_n}\]

Parseval’s Theorem is a special case of the power theorem:

\[\sum_n |x_n|^2 = \frac{1}{N} \sum_n |\tilde x_n|^2\]

When the fourier space array contains only half the spectrum, making use of hermitian symmetry, extra care has to be taken when performing the sum.

# TODO

muFFT fourier transform:

fft and ifft never applies the normalisation factor, meaning that you will need to multiply ifft(fft) by 1 / np.prod(nb_grid_pts) = fftengine.normalisation) in order to have a roundtrip.

muFFT vs. np.fft:

Normalisation:

np.fft.rfft <–> fftengine.fft

np.fft.irfft <–> fftengine.ifft * fftengine.normalisation

2D FFT:

numpy by default transforms the last index first.

muFFT the first ` real_buffer.array()[..] = a fftengine.fft(real_buffer, fourier_buffer) fourier_buffer <--> np.rfft2(a.T).T <--> np.fft.rfft2(a, axes=(1,0)) ` # FIXME: @pastewka: I expected the fourier array to be transposed, so there is a # wrapper swapping the indexes and the array # is transposed in memory ?

class ContactMechanics.FFTElasticHalfSpace.BndSet(large, small)

Bases: tuple

Methods

`count`(value, /)	Return number of occurrences of value.
`index`(value[, start, stop])	Return first index of value.

large: Alias for field number 0

small: Alias for field number 1

class ContactMechanics.FFTElasticHalfSpace.FreeFFTElasticHalfSpace(**kwargs)

Bases: PeriodicFFTElasticHalfSpace

Uses the FFT to solve the displacements and stresses in an non-periodic elastic Halfspace due to a given array of point forces. Uses the Green’s functions formulaiton of Johnson (1985, p. 54). The application of the FFT to a nonperiodic domain is explained in Hockney (1969, p. 178.)

K. L. Johnson. (1985). Contact Mechanics. [Online]. Cambridge: Cambridge University Press. Available from: Cambridge Books Online <http://dx.doi.org/10.1017/CBO9781139171731> [Accessed 16 February 2015]

R. W. HOCKNEY, “The potential calculation and some applications,” Methods of Computational Physics, B. Adler, S. Fernback and M. Rotenberg (Eds.), Academic Press, New York, 1969, pp. 136-211.

Attributes:

area_per_pt
communicator: Return the MPI communicator
dim: return the substrate’s physical dimension
domain_boundary_mask: Returns a mask of the points that are on the boundary of the domain
fourier_locations: When working in Parallel one processor holds only Part of the Data
fourier_slices: When working in Parallel one processor holds only Part of the Data
local_topography_subdomain_slices: slice representing the local subdomain without the padding area
nb_domain_grid_pts: usually, the nb_grid_pts of the system is equal to the geometric
nb_fourier_grid_pts: When working in Parallel one processor holds only Part of the Data
nb_grid_pts
nb_subdomain_grid_pts: When working in Parallel one processor holds only Part of the Data
physical_sizes
subdomain_locations: When working in Parallel one processor holds only Part of the Data
subdomain_slices: When working in Parallel one processor holds only Part of the Data
topography_nb_subdomain_grid_pts
topography_subdomain_locations
topography_subdomain_slices

Methods

`FreeBoundaryError`(message)	called when the forces overlap into the padding region (i.e. the outer ring of the force array equals zero), needing an increase of the nb_grid_pts.
`check`([force])	Checks wether force is still in the value range handled correctly :param force:
`check_boundaries`([force, tol])	Raises an error if the forces are not zero at the boundary of the active domain
`compute`(disp[, pot, forces])	computes and stores the elastic energy and/or surface forces the as function of the surface displacement.
`compute_k`(disp_k[, pot, forces])	computes and stores the elastic energy and/or surface forces the as function of the surface displacement in Fourier Space.
`evaluate`(disp[, pot, forces])	Evaluates the elastic energy and the point forces Keyword Arguments: disp -- array of distances pot -- (default True) if true, returns potential energy forces -- (default False) if true, returns forces
`evaluate_disp`(forces)	Computes the displacement due to a given force array Keyword Arguments: forces -- a numpy array containing point forces (not pressures)
`evaluate_elastic_energy`(forces, disp)	computes and returns the elastic energy due to forces and displacements Arguments: forces -- array of forces disp -- array of displacements
`evaluate_elastic_energy_k_space`(kforces, kdisp)	Computes the Energy due to forces and displacements using their Fourier representation.
`evaluate_force`(disp)	Computes the force (not pressures) due to a given displacement array.
`evaluate_k`(disp_k[, pot, forces])	Evaluates the elastic energy and the point forces in fourier sapce or k-space or reciprocal space.
`evaluate_k_disp`(forces)	Computes the K-space displacement due to a given force array
`evaluate_k_force`(disp)	Computes the K-space forces (not pressures) due to a given displacement array.
`evaluate_k_force_k`(disp_k)	Computes the K-space forces (not pressures) due to a given K-space displacement array.
`evaluate_scalar_product_k_space`(ka, kb)	Computes the scalar product, i.e. the power, between the a and b, given their fourier representation.
`is_periodic`()	non-periodic substrates can use some optimisations
`spawn_child`(nb_grid_pts)	returns an instance with same physical properties with a smaller computational grid

Error

exception FreeBoundaryError(message)

Bases: Exception

called when the forces overlap into the padding region (i.e. the outer ring of the force array equals zero), needing an increase of the nb_grid_pts

check(force=None): Checks wether force is still in the value range handled correctly :param force:

check_boundaries(force=None, tol=0)

Raises an error if the forces are not zero at the boundary of the active domain

Parameters:: force

property domain_boundary_mask

Returns a mask of the points that are on the boundary of the domain

Return type:: bool ndarray with size self.topography_nb_subdomain_grid_pts

evaluate_disp(forces)

Computes the displacement due to a given force array Keyword Arguments: forces – a numpy array containing point forces (not pressures)

if running in MPI this should be only the forces in the Subdomain

if running in serial one can give the force array with or without the padded region

name = 'free_fft_elastic_halfspace'

property nb_domain_grid_pts: usually, the nb_grid_pts of the system is equal to the geometric nb_grid_pts (of the surface). For example free boundary conditions, require the computational nb_grid_pts to differ from the geometric one, see FreeFFTElasticHalfSpace.

spawn_child(nb_grid_pts): returns an instance with same physical properties with a smaller computational grid

property topography_nb_subdomain_grid_pts

class ContactMechanics.FFTElasticHalfSpace.PeriodicFFTElasticHalfSpace(**kwargs)

Bases: ElasticSubstrate

Uses the FFT to solve the displacements and stresses in an elastic Halfspace due to a given array of point forces. This halfspace implementation cheats somewhat: since a net pressure would result in infinite displacement, the first term of the FFT is systematically dropped. The implementation follows the description in Stanley & Kato J. Tribol. 119(3), 481-485 (Jul 01, 1997)

Attributes:

area_per_pt
communicator: Return the MPI communicator
dim: return the substrate’s physical dimension
fourier_locations: When working in Parallel one processor holds only Part of the Data
fourier_slices: When working in Parallel one processor holds only Part of the Data
local_topography_subdomain_slices: slice representing the local subdomain without the padding area
nb_domain_grid_pts: usually, the nb_grid_pts of the system is equal to the geometric
nb_fourier_grid_pts: When working in Parallel one processor holds only Part of the Data
nb_grid_pts
nb_subdomain_grid_pts: When working in Parallel one processor holds only Part of the Data
physical_sizes
subdomain_locations: When working in Parallel one processor holds only Part of the Data
subdomain_slices: When working in Parallel one processor holds only Part of the Data
topography_nb_subdomain_grid_pts
topography_subdomain_locations
topography_subdomain_slices

Methods

`check`([force])	Checks wether force is still in the value range handled correctly.
`compute`(disp[, pot, forces])	computes and stores the elastic energy and/or surface forces the as function of the surface displacement.
`compute_k`(disp_k[, pot, forces])	computes and stores the elastic energy and/or surface forces the as function of the surface displacement in Fourier Space.
`evaluate`(disp[, pot, forces])	Evaluates the elastic energy and the point forces Keyword Arguments: disp -- array of distances pot -- (default True) if true, returns potential energy forces -- (default False) if true, returns forces
`evaluate_disp`(forces)	Computes the displacement due to a given force array Keyword Arguments: forces -- a numpy array containing point forces (not pressures)
`evaluate_elastic_energy`(forces, disp)	computes and returns the elastic energy due to forces and displacements Arguments: forces -- array of forces disp -- array of displacements
`evaluate_elastic_energy_k_space`(kforces, kdisp)	Computes the Energy due to forces and displacements using their Fourier representation.
`evaluate_force`(disp)	Computes the force (not pressures) due to a given displacement array.
`evaluate_k`(disp_k[, pot, forces])	Evaluates the elastic energy and the point forces in fourier sapce or k-space or reciprocal space.
`evaluate_k_disp`(forces)	Computes the K-space displacement due to a given force array
`evaluate_k_force`(disp)	Computes the K-space forces (not pressures) due to a given displacement array.
`evaluate_k_force_k`(disp_k)	Computes the K-space forces (not pressures) due to a given K-space displacement array.
`evaluate_scalar_product_k_space`(ka, kb)	Computes the scalar product, i.e. the power, between the a and b, given their fourier representation.
`is_periodic`()	non-periodic substrates can use some optimisations
`spawn_child`(dummy)	does nothing for most substrates.

Error

property area_per_pt

property communicator: Return the MPI communicator

property dim: return the substrate’s physical dimension

evaluate(disp, pot=True, forces=False): Evaluates the elastic energy and the point forces Keyword Arguments: disp – array of distances pot – (default True) if true, returns potential energy forces – (default False) if true, returns forces

evaluate_disp(forces): Computes the displacement due to a given force array Keyword Arguments: forces – a numpy array containing point forces (not pressures)

evaluate_elastic_energy(forces, disp): computes and returns the elastic energy due to forces and displacements Arguments: forces – array of forces disp – array of displacements

evaluate_elastic_energy_k_space(kforces, kdisp)

Computes the Energy due to forces and displacements using their Fourier representation.

\[ \begin{align}\begin{aligned}E_{el} &= - \frac{1}{2} \sum_{ij} u_{ij} f_{ij}\\ &= - \frac{1}{2} \frac{1}{n_x n_y} \sum_{kl} \tilde u_{kl} \overline{\tilde f_{kl}}\end{aligned}\end{align} \]

(\(\tilde f_{ij} = - \tilde K_{ijkl} \tilde u_{kl}\))

In a parallelized code kforces and kdisp contain only the slice attributed to this processor

Parameters:

kforces – array of complex type and of size substrate.nb_fourier_grid_pts Fourier representation (output of a 2D rfftn) of the forces acting on the grid points
kdisp – array of complex type and of physical_sizes substrate.nb_fourier_grid_pts Fourier representation (output of a 2D rfftn) of the displacements of the grid points

Returns:

The elastic energy due to the forces and displacements

Return type:

E

evaluate_force(disp)

Computes the force (not pressures) due to a given displacement array.

Keyword Arguments: disp – a numpy array containing point displacements

evaluate_k(disp_k, pot=True, forces=False): Evaluates the elastic energy and the point forces in fourier sapce or k-space or reciprocal space.

evaluate_k_disp(forces)

Computes the K-space displacement due to a given force array

Parameters:: forces (ndarray) – a numpy array containing point forces (not pressures)
Returns:: displacement – displacement in k-space
Return type:: ndarray

evaluate_k_force(disp)

Computes the K-space forces (not pressures) due to a given displacement array.

Keyword Arguments: disp – a numpy array containing point displacements

evaluate_k_force_k(disp_k)

Computes the K-space forces (not pressures) due to a given K-space displacement array.

Parameters:: disp (ndarray k-space) – a numpy k-space array containing point displacements
Returns:: force_k
Return type:: nd array k-sapce forces

evaluate_scalar_product_k_space(ka, kb)

Computes the scalar product, i.e. the power, between the a and b, given their fourier representation.

Power theorem:

\[P = \sum_{ij} a_{ij} b_{ij} = \frac{1}{n_x n_y}\sum_{ij} \tilde a_{ij} \overline{\tilde b_{ij}}\]

Note that for a, b real,

\[P = \sum_{kl} Re(\tilde a_{kl}) Re(\tilde b_{kl})\]

Im(tilde a_{kl}) Im(tilde b_{kl})

Parameters:

ka – arrays of complex type and of size substrate.nb_fourier_grid_pts Fourier representation (output of a 2D rfftn) a (resp. b) (nx, ny real array)
kb – arrays of complex type and of size substrate.nb_fourier_grid_pts Fourier representation (output of a 2D rfftn) a (resp. b) (nx, ny real array)

Returns:

The scalar product of a and b

Return type:

P

property fourier_locations

When working in Parallel one processor holds only Part of the Data

Returns:

property fourier_slices

When working in Parallel one processor holds only Part of the Data

Returns:

property local_topography_subdomain_slices: slice representing the local subdomain without the padding area

name = 'periodic_fft_elastic_halfspace'

property nb_domain_grid_pts: usually, the nb_grid_pts of the system is equal to the geometric nb_grid_pts (of the surface). For example free boundary conditions, require the computational nb_grid_pts to differ from the geometric one, see FreeFFTElasticHalfSpace.

property nb_fourier_grid_pts

When working in Parallel one processor holds only Part of the Data

Returns:

property nb_grid_pts

property nb_subdomain_grid_pts

When working in Parallel one processor holds only Part of the Data

Returns:

property physical_sizes

property subdomain_locations

When working in Parallel one processor holds only Part of the Data

Returns:

property subdomain_slices

When working in Parallel one processor holds only Part of the Data

Returns:

property topography_nb_subdomain_grid_pts

property topography_subdomain_locations

property topography_subdomain_slices

ContactMechanics.Factory module

Implements a convenient Factory function for Contact System creation

ContactMechanics.Factory.make_contact_system(topography, *args, **kwargs)

Factory function for contact systems. Checks the compatibility between the substrate, interaction method and surface and returns an object of the appropriate type to handle it. The returned object is always of a subtype of SystemBase.

This function can be called in different ways: - provide the elastic propoerties of the substrates

>>> topography.make_contact_system(young=1)
for other possibilities, i.e. finite thickness\
check the arguments of `FFTElasticHalfspace` or `FreeFFTElasticHalfspace`

provide an already instantiated substrate >>> topography.make_contact_system(substrate=substrate)

ContactMechanics.Factory.make_plastic_contact_system(topography, *args, **kwargs)

ContactMechanics.Factory.make_plastic_system(*args, **kwargs): Factory function for contact systems. Checks the compatibility between the substrate, interaction method and surface and returns an object of the appropriate type to handle it. The returned object is always of a subtype of SystemBase.

ContactMechanics.Factory.make_system(*args, **kwargs): Factory function for contact systems. Checks the compatibility between the substrate, interaction method and surface and returns an object of the appropriate type to handle it. The returned object is always of a subtype of SystemBase.

ContactMechanics.PipelineFunction module

ContactMechanics.PipelineFunction.contact_mechanics(self, substrate=None, nsteps=None, offsets=None, pressures=None, hardness=None, maxiter=100, results_callback=None, optimizer_kwargs={})

Carry out an automated contact mechanics calculations. The pipeline function return thermodynamic data (averages over the contact area, e.g. the total force or the total area). Spatially resolved data (pressure maps, displacement maps, etc.) are passed to the callback function. If this data is reqired, the callback function needs to take care of analyzing or storing it.

Parameters:

self (SurfaceTopography.UniformTopographyInterface) – Topography on which to carry out the contact calculation.
substrate (str, optional) – Specifies whether substrate should be ‘periodic’ or ‘nonperiodic’. If set to None, it will be chosen according to whether the topography is periodic or nonperiodic. (Default: None)
nsteps (int, optional) – Number of contact steps. (Default: 10)
offsets (list of floats, optional) – List with offsets. Can only be set if nsteps and pressures is set to None. (Default: None)
pressures (list of floats, optional) – List with pressures in units of E*. Can only be set if nsteps and offsets is set to None. (Default: None)
hardness (float, optional) – Hardness in units of E*. Calculation is fully elastic if set to None. (Default: None)
maxiter (int, optional) – Maximum number of interations. (Default: 100)
results_callback (func, optional) – Callback function receiving displacement, pressure, etc. fields. (Default: None)
optimizer_kwargs (dict, optional) – Optional arguments passed on to the optimizer. (Default: {})

Returns:

mean_pressure (np.ndarray) – Array with mean pressure for each calculation step.
total_contact_area (np.ndarray) – Array with total area for each calculation step.
mean_displacement (np.ndarray) – Array with mean displacement for each calculation step.
mean_gap (np.ndarray) – Array with mean gap for each calculation step.
converged (np.ndarray) – Convergence information for each calculation step. Unconverged results are still returned but should be interpreted with care.

ContactMechanics.PlasticSystemSpecialisations module

implements plastic mapping algorithms for contact systems

class ContactMechanics.PlasticSystemSpecialisations.PlasticNonSmoothContactSystem(substrate, surface)

Bases: NonSmoothContactSystem

This system implements a simple penetration hardness model.

Attributes:

nb_grid_pts: For systems, nb_grid_pts can become non-trivial

Methods

`callback`([force])	Simple callback function that can be handed over to scipy's minimize to get updates during minimisation Parameters: force -- (default False) whether to include the norm of the force vector in the update message
`compute_contact_area`()	computes and returns the total contact area
`compute_contact_coordinates`()	returns an array of all coordinates, where contact pressure is repulsive.
`compute_gap`(disp, offset, *profile_args, ...)	evaluate the gap between surface and substrate.
`compute_nb_contact_pts`()	compute and return the number of contact points.
`compute_normal_force`()	computes and returns the sum of all forces
`compute_relative_contact_area`()	compute and return the relative contact area:
`dual_hessian_product`(pressure)	Returns the hessian product of the dual_objective function.
`dual_minimize_proxy`(offset[, init_force, ...])	Convenience function for DUAL minimisation (pixel forces as variables).
`dual_objective`(offset[, gradient])	Objective function to handle dual objective, i.e. the Legendre transformation from displacements as variable to pressures (the Lagrange multiplier) as variable.
`evaluate`(disp, offset[, pot, forces, logger])	Compute the energies and forces in the system for a given displacement field.
`evaluate_dual`(press, offset[, forces])	Computes the energies and forces in the system for a given pressure field.
`handles`(substrate_type, surface_type, ...)	determines whether this class can handle the proposed system composition Keyword Arguments: substrate_type -- instance of ElasticSubstrate subclass surface_type --
`hessian_product`(disp)	Computes the Hessian product for the objective function.
`is_proxy`()	subclasses may not be able to implement the full interface because they try to do something smart and internally compute a different system. They should declare to to be proxies and provide a method called cls. deproxyfied() that returns the energy, force and displacement of the full problem based on its internal state. E.g at the end of an optimization, you could have: if system.is_proxy(): energy, force, disp = system.deproxyfied().
`logger_input`()	Describes the current state of the system (during minimization)
`objective`(offset[, disp0, gradient, logger])	This helper method exposes a scipy.optimize-friendly interface to the evaluate() method.
`primal_hessian_product`(gap)	Returns the hessian product of the primal_objective function.
`primal_minimize_proxy`(offset[, init_gap, ...])	This function is a convenience function that simplifies the process of solving the primal minimisation problem where the gap is the variable.
`primal_objective`(offset[, gradient])	Solves the primal objective using gap as the variable.
`shape_minimisation_input`(in_array)	For minimisation of smart systems, the initial guess array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues).
`shape_minimisation_output`(in_array)	For minimisation of smart systems, the output array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues).

minimize_proxy

static handles(substrate_type, surface_type, is_domain_decomposed): determines whether this class can handle the proposed system composition Keyword Arguments: substrate_type – instance of ElasticSubstrate subclass surface_type –

minimize_proxy(**kwargs)

ContactMechanics.Substrates module

Base class for continuum mechanics models of halfspaces

class ContactMechanics.Substrates.ElasticSubstrate

Bases: Substrate

Generic baseclass for elastic substrates

Attributes:

communicator: Return the MPI communicator
nb_domain_grid_pts
nb_subdomain_grid_pts: When working in Parallel one processor holds only Part of the Data
subdomain_locations: When working in Parallel one processor holds only Part of the Data

Methods

`check`([force])	Checks wether force is still in the value range handled correctly.
`compute`(disp[, pot, forces])	computes and stores the elastic energy and/or surface forces the as function of the surface displacement.
`compute_k`(disp_k[, pot, forces])	computes and stores the elastic energy and/or surface forces the as function of the surface displacement in Fourier Space.
`evaluate`(disp[, pot, forces])	computes and returns the elastic energy and/or surface forces as function of the surface displacement.
`is_periodic`()	non-periodic substrates can use some optimisations
`spawn_child`(dummy)	does nothing for most substrates.

Error

compute(disp, pot=True, forces=False): computes and stores the elastic energy and/or surface forces the as function of the surface displacement. Note that forces, not surface pressures are expected. This is contrary to most formulations in the literature, but convenient in the code (consistency with the softWall interaction potentials). This choice may come back to bite me. Parameters: gap – array containing the point-wise gap values pot – (default True) whether the energy should be evaluated forces – (default False) whether the forces should be evaluated

compute_k(disp_k, pot=True, forces=False): computes and stores the elastic energy and/or surface forces the as function of the surface displacement in Fourier Space. Note that forces, not surface pressures are expected. This is contrary to most formulations in the literature, but convenient in the code (consistency with the softWall interaction potentials). This choice may come back to bite me. Parameters: gap – array containing the point-wise gap values pot – (default True) whether the energy should be evaluated forces – (default False) whether the forces should be evaluated

evaluate(disp, pot=True, forces=False): computes and returns the elastic energy and/or surface forces as function of the surface displacement. See docstring for ‘compute’ for more details Parameters: gap – array containing the point-wise gap values pot – (default True) whether the energy should be evaluated forces – (default False) whether the forces should be evaluated

name = 'generic_elastic_halfspace'

class ContactMechanics.Substrates.PlasticSubstrate

Bases: Substrate

Generic baseclass for plastic substrates

Attributes:

communicator: Return the MPI communicator
nb_domain_grid_pts
nb_subdomain_grid_pts: When working in Parallel one processor holds only Part of the Data
subdomain_locations: When working in Parallel one processor holds only Part of the Data

Methods

`check`([force])	Checks wether force is still in the value range handled correctly.
`is_periodic`()	non-periodic substrates can use some optimisations
`spawn_child`(dummy)	does nothing for most substrates.

Error

name = 'generic_plastic_halfspace'

class ContactMechanics.Substrates.Substrate

Bases: object

Generic baseclass from which all substate classes derive

Attributes:

communicator: Return the MPI communicator
nb_domain_grid_pts
nb_subdomain_grid_pts: When working in Parallel one processor holds only Part of the Data
subdomain_locations: When working in Parallel one processor holds only Part of the Data

Methods

`check`([force])	Checks wether force is still in the value range handled correctly.
`is_periodic`()	non-periodic substrates can use some optimisations
`spawn_child`(dummy)	does nothing for most substrates.

Error

exception Error: Bases: Exception

check(force=None): Checks wether force is still in the value range handled correctly. In this case all forces are ok. :param force:

abstract property communicator: Return the MPI communicator

classmethod is_periodic(): non-periodic substrates can use some optimisations

name = 'generic_halfspace'

abstract property nb_domain_grid_pts

abstract property nb_subdomain_grid_pts

When working in Parallel one processor holds only Part of the Data

Returns:

spawn_child(dummy): does nothing for most substrates.

abstract property subdomain_locations

When working in Parallel one processor holds only Part of the Data

Returns:

ContactMechanics.Systems module

Defines the interface for contact mechanics systems

exception ContactMechanics.Systems.IncompatibleFormulationError: Bases: Exception

exception ContactMechanics.Systems.IncompatibleResolutionError: Bases: Exception

class ContactMechanics.Systems.NonSmoothContactSystem(substrate, surface)

Bases: SystemBase

For non-smooth contact mechanics (i.e, the equlibrium is the solution to a constrained optimisation problem with a non-zero gradient of the energy functional at the solution). The classic contact problems, for which the interaction between the tribopartners is just non-penetration without adhesion, belong to this type of system

Attributes:

nb_grid_pts: For systems, nb_grid_pts can become non-trivial

Methods

`callback`([force])	Simple callback function that can be handed over to scipy's minimize to get updates during minimisation Parameters: force -- (default False) whether to include the norm of the force vector in the update message
`compute_contact_area`()	computes and returns the total contact area
`compute_contact_coordinates`()	returns an array of all coordinates, where contact pressure is repulsive.
`compute_gap`(disp, offset, *profile_args, ...)	evaluate the gap between surface and substrate.
`compute_nb_contact_pts`()	compute and return the number of contact points.
`compute_normal_force`()	computes and returns the sum of all forces
`compute_relative_contact_area`()	compute and return the relative contact area:
`dual_hessian_product`(pressure)	Returns the hessian product of the dual_objective function.
`dual_minimize_proxy`(offset[, init_force, ...])	Convenience function for DUAL minimisation (pixel forces as variables).
`dual_objective`(offset[, gradient])	Objective function to handle dual objective, i.e. the Legendre transformation from displacements as variable to pressures (the Lagrange multiplier) as variable.
`evaluate`(disp, offset[, pot, forces, logger])	Compute the energies and forces in the system for a given displacement field.
`evaluate_dual`(press, offset[, forces])	Computes the energies and forces in the system for a given pressure field.
`handles`(substrate_type, surface_type, ...)	determines whether this class can handle the proposed system composition Keyword Arguments: substrate_type -- instance of ElasticSubstrate subclass surface_type --
`hessian_product`(disp)	Computes the Hessian product for the objective function.
`is_proxy`()	subclasses may not be able to implement the full interface because they try to do something smart and internally compute a different system. They should declare to to be proxies and provide a method called cls. deproxyfied() that returns the energy, force and displacement of the full problem based on its internal state. E.g at the end of an optimization, you could have: if system.is_proxy(): energy, force, disp = system.deproxyfied().
`logger_input`()	Describes the current state of the system (during minimization)
`minimize_proxy`([solver])	Convenience function.
`objective`(offset[, disp0, gradient, logger])	This helper method exposes a scipy.optimize-friendly interface to the evaluate() method.
`primal_hessian_product`(gap)	Returns the hessian product of the primal_objective function.
`primal_minimize_proxy`(offset[, init_gap, ...])	This function is a convenience function that simplifies the process of solving the primal minimisation problem where the gap is the variable.
`primal_objective`(offset[, gradient])	Solves the primal objective using gap as the variable.
`shape_minimisation_input`(in_array)	For minimisation of smart systems, the initial guess array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues).
`shape_minimisation_output`(in_array)	For minimisation of smart systems, the output array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues).

compute_contact_coordinates(): returns an array of all coordinates, where contact pressure is repulsive. Useful for evaluating the number of contact islands etc.

compute_nb_contact_pts(): compute and return the number of contact points. Note that this is of no physical interest, as it is a purely numerical artefact

compute_normal_force(): computes and returns the sum of all forces

dual_hessian_product(pressure)

Returns the hessian product of the dual_objective function.

The hessian product is the result of applying the hessian matrix (second derivatives of the objective function) to the pressure vector. This is used in optimization algorithms that utilize second-order information.

Parameters:: pressure (array_like) – The pressure vector to which the hessian matrix is applied.
Returns:: hessp – The hessian product, which is the result of applying the hessian matrix to the pressure vector.
Return type:: array_like

dual_minimize_proxy(offset, init_force=None, solver='ccg-without-restart', gtol=1e-08, maxiter=1000)

Convenience function for DUAL minimisation (pixel forces as variables). This function simplifies the process of solving the dual minimisation problem by encapsulating the use of constrained minimisation.

Parameters:

offset (float) – Determines the indentation depth.
init_force (array_like, optional) – Initial guess for the force. If not provided, it defaults to None.
solver (str, optional) – The solver to be used for the minimisation. It can be one of ‘ccg-without-restart’, ‘ccg-with-restart’, or ‘l-bfgs-b’. If not provided, it defaults to ‘ccg-without-restart’.
gtol (float, optional) – The gradient tolerance for the solver. If not provided, it defaults to 1e-8.
maxiter (int, optional) – The maximum number of iterations allowed for the solver to converge. If not provided, it defaults to 1000.

Returns:

result – The result of the minimisation. It contains information about the optimisation result, including the final gap, force, and displacement of the system at the solution.

Return type:

OptimizeResult

dual_objective(offset, gradient=True)

Objective function to handle dual objective, i.e. the Legendre transformation from displacements as variable to pressures (the Lagrange multiplier) as variable.

Parameters:

offset (float) – Constant value to add to the surface heights.
gradient (bool, optional) – Whether to return the gradient in addition to the energy. Default is True.

Returns:

energy (float) – Value of total energy.
gradient (array_like) – Value of the gradient (array) or the value of gap (if gradient is True).

Notes

The objective function is defined as:

\[\begin{split}\min_\lambda \ q(\lambda) = \frac{1}{2}\lambda_i K^{-1}_{ij} \lambda_j - \lambda_i h_i \\ \\ \nabla q = K^{-1}_{ij} \lambda_j - h_i \hspace{0.1cm} \text{which is,} \\ \text{gap} = \text{displacement} - \text{height} \\\end{split}\]

evaluate(disp, offset, pot=True, forces=False, logger=None)

Compute the energies and forces in the system for a given displacement field.

This method calculates the gap between the surface and the substrate by calling the compute_gap method. It then computes the displacement field by calling the compute method of the substrate.

If potential energy is to be computed, it is set to the energy of the substrate. Otherwise, it is set to None.

If forces are to be computed, they are set to the force of the substrate. Otherwise, they are set to None.

If a logger is provided, it logs the current state of the system.

Parameters:

disp (array_like) – The displacement field for which the energies and forces are to be computed.
offset (float) – The offset value to be used in the computation.
pot (bool, optional) – If True, the potential energy in the system is also computed. Default is True.
forces (bool, optional) – If True, the forces in the system are also computed. Default is False.
logger (Logger, optional) – Logger object to log information at each iteration. Default is None.

Returns:

energy (float) – Total energy of the system. If potential energy is not computed, it is None.
force (array_like) – Forces in the system. If forces are not computed, they are None.

evaluate_dual(press, offset, forces=False)

Computes the energies and forces in the system for a given pressure field.

This method calculates the displacement field corresponding to the given pressure field by calling the evaluate_disp method of the substrate. The negative of the pressure field is passed as an argument to the evaluate_disp method.

If forces are to be computed, the gradient is calculated as the difference between the displacement and the sum of the surface heights and the offset. Otherwise, the gradient is set to None.

The energy is then computed as half the sum of the product of the pressure and displacement fields, minus the sum of the product of the pressure and the sum of the surface heights and the offset.

Parameters:

press (array_like) – The pressure field for which the displacement field is to be computed.
offset (float) – The offset value to be used in the computation.
forces (bool, optional) – If True, the forces in the system are also computed. Default is False.

Returns:

energy (float) – Total energy of the system.
gradient (array_like) – Gradient, which is the difference between the displacement and the sum of the surface heights and the offset. If forces are not computed, the gradient is None.

static handles(substrate_type, surface_type, is_domain_decomposed): determines whether this class can handle the proposed system composition Keyword Arguments: substrate_type – instance of ElasticSubstrate subclass surface_type –

hessian_product(disp)

Computes the Hessian product for the objective function.

This method calculates the Hessian product by calling the primal_hessian_product method. The Hessian product is the result of applying the Hessian matrix (second derivatives of the objective function) to the displacement vector. This is used in optimization algorithms that utilize second-order information.

Parameters:: disp (array_like) – The displacement vector to which the Hessian matrix is applied. It can be an array of shape nb_subdomain_grid_pts or a flattened version of it.
Returns:: The Hessian product, which is the result of applying the Hessian matrix to the displacement vector.
Return type:: array_like

logger_input()

Describes the current state of the system (during minimization)

Output is suited to be passed to ContactMechanics.Tools.Logger.Logger

Returns:

headers (list of strings)
values (list)

minimize_proxy(solver=<function doi.__call__.<locals>.func_with_doi>, **kwargs): Convenience function. Eliminates boilerplate code for most minimisation problems by encapsulating the use of constrained minimisation.

property nb_grid_pts: For systems, nb_grid_pts can become non-trivial

objective(offset, disp0=None, gradient=False, logger=None)

This helper method exposes a scipy.optimize-friendly interface to the evaluate() method. It is used for optimization purposes and ensures that the shape of displacement is maintained. It also allows setting the offset and ‘forces’ flag without using scipy’s argument passing interface.

Parameters:

offset (float) – Determines the indentation depth.
disp0 (array_like, optional) – Unused variable, present only for interface compatibility with inheriting classes. Default is None.
gradient (bool, optional) – If True, the gradient is supposed to be used. Default is False.
logger (Logger, optional) – Logger object to log information at each iteration. Default is None.

Returns:

fun – A function of only displacement. If gradient is True, this function returns the energy and the negative of the force when called with displacement. If gradient is False, it returns only the energy.

Return type:

function

primal_hessian_product(gap)

Returns the hessian product of the primal_objective function.

The hessian product is the result of applying the hessian matrix (second derivatives of the objective function) to the gap vector. This is used in optimization algorithms that utilize second-order information.

Parameters:: gap (array_like) – The gap vector to which the hessian matrix is applied.
Returns:: hessp – The hessian product, which is the result of applying the hessian matrix to the gap vector.
Return type:: array_like

Notes

The hessian product is computed as the negative of the force evaluated at the reshaped gap vector.

primal_minimize_proxy(offset, init_gap=None, solver='ccg-without-restart', gtol=1e-08, maxiter=1000)

This function is a convenience function that simplifies the process of solving the primal minimisation problem where the gap is the variable. It does this by encapsulating the use of constrained minimisation.

Parameters:

offset (float) – This parameter determines the indentation depth.
init_gap (array_like, optional) – This is the initial guess for the gap. If not provided, it defaults to None.
solver (str, optional) – This is the solver to be used for the minimisation. It can be one of ‘ccg-without-restart’, ‘ccg-with-restart’, or ‘l-bfgs-b’. If not provided, it defaults to ‘ccg-without-restart’.
gtol (float, optional) – This is the gradient tolerance for the solver. If not provided, it defaults to 1e-8.
maxiter (int, optional) – This is the maximum number of iterations allowed for the solver to converge. If not provided, it defaults to 1000.

Returns:

result – The result of the minimisation. It contains information about the optimisation result, including the final gap, force, and displacement of the system at the solution.

Return type:

OptimizeResult

primal_objective(offset, gradient=True)

Solves the primal objective using gap as the variable. This function can be fed directly to standard solvers such as scipy solvers etc. and returns the elastic energy and its gradient (negative of the forces) as a function of the gap.

Parameters:

offset (float) – Constant value to add to the surface heights.
gradient (bool, optional) – Return gradient in addition to the energy. Default is True.

Returns:

energy (float) – Value of total energy.
force (array_like) – Value of the forces per surface node (only if gradient is True).

Notes

The objective function is defined as:

\[\begin{split}\min_u f(u) = 1/2 u_i K_{ij} u_j \\ \\ \nabla f = K_{ij} u_j \ \ \ \text{which is the Force.} \\\end{split}\]

class ContactMechanics.Systems.SystemBase(substrate, surface)

Bases: object

Base class for contact systems

Attributes:

nb_grid_pts: For systems, nb_grid_pts can become non-trivial

Methods

`callback`([force])	Simple callback function that can be handed over to scipy's minimize to get updates during minimisation Parameters: force -- (default False) whether to include the norm of the force vector in the update message
`compute_contact_area`()	computes and returns the total contact area
`compute_gap`(disp, offset, *profile_args, ...)	evaluate the gap between surface and substrate.
`compute_nb_contact_pts`()	compute and return the number of contact points.
`compute_normal_force`()	evaluates and returns the normal force between substrate and surface
`compute_relative_contact_area`()	compute and return the relative contact area:
`evaluate`(disp, offset[, pot, forces])	Compute the energies and forces in the system for a given displacement field
`handles`(substrate_type, surface_type, ...)	returns whether this class (in practice a subclass) handles this combination of types Keyword Arguments: substrate_type -- self-explanatory surface_type -- self-explanatory is_domain_decomposed: some systems cannot handle parallel computation
`is_proxy`()	subclasses may not be able to implement the full interface because they try to do something smart and internally compute a different system. They should declare to to be proxies and provide a method called cls. deproxyfied() that returns the energy, force and displacement of the full problem based on its internal state. E.g at the end of an optimization, you could have: if system.is_proxy(): energy, force, disp = system.deproxyfied().
`minimize_proxy`([offset, ...])	Convenience function.
`objective`(offset[, disp0, gradient, logger])	This helper method exposes a scipy.optimize-friendly interface to the evaluate() method. Use this for optimization purposes, it makes sure that the shape of disp is maintained and lets you set the offset and 'forces' flag without using scipy's cumbersome argument passing interface. Returns a function of only disp Keyword Arguments: offset -- determines indentation depth disp0 -- preexisting displacement. influences e.g., the physical_sizes of the proxy system in some 'smart' system subclasses gradient -- (default False) whether the gradient is supposed to be used logger -- (default None) log information at every iteration.
`shape_minimisation_input`(in_array)	For minimisation of smart systems, the initial guess array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues).
`shape_minimisation_output`(in_array)	For minimisation of smart systems, the output array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues).

callback(force=False): Simple callback function that can be handed over to scipy’s minimize to get updates during minimisation Parameters: force – (default False) whether to include the norm of the force

vector in the update message

compute_contact_area(): computes and returns the total contact area

compute_gap(disp, offset, *profile_args, **profile_kwargs): evaluate the gap between surface and substrate. Convention is that non-penetrating contact has gap >= 0

abstract compute_nb_contact_pts(): compute and return the number of contact points. Note that this is of no physical interest, as it is a purely numerical artefact

abstract compute_normal_force(): evaluates and returns the normal force between substrate and surface

compute_relative_contact_area()

compute and return the relative contact area:: A
Aᵣ = ──: A₀

abstract evaluate(disp, offset, pot=True, forces=False): Compute the energies and forces in the system for a given displacement field

static handles(substrate_type, surface_type, is_domain_decomposed): returns whether this class (in practice a subclass) handles this combination of types Keyword Arguments: substrate_type – self-explanatory surface_type – self-explanatory is_domain_decomposed: some systems cannot handle parallel computation

classmethod is_proxy(): subclasses may not be able to implement the full interface because they try to do something smart and internally compute a different system. They should declare to to be proxies and provide a method called cls. deproxyfied() that returns the energy, force and displacement of the full problem based on its internal state. E.g at the end of an optimization, you could have: if system.is_proxy():

energy, force, disp = system.deproxyfied()

minimize_proxy(offset=0, initial_displacements=None, method='L-BFGS-B', gradient=True, lbounds=None, ubounds=None, callback=None, logger=None, **kwargs)

Convenience function. Eliminates boilerplate code for most minimisation problems by encapsulating the use of scipy.minimize for common default options. In the case of smart proxy systems, this may also encapsulate things like dynamics computation of safety margins, extrapolation of results onto the proxied system, etc.

Parameters: offset : float

determines indentation depth

initial_displacements(default zero)

initial guess for displacement field. If not chosen appropriately, results may be unreliable.

methodstring or callable

(defaults to L-BFGS-B, see scipy documentation). Be sure to choose method that can handle high-dimensional parameter spaces.

optionsdict

(default None) options to be passed to the minimizer method

gradientbool

(default True) whether to use the gradient or not

lboundsarray of shape substrate.subdomain_nb_grid_pts or

substrate.topography_subdomain_nb_grid_pts or “auto”

(default None): nodal ceiling/floor

uboundsarray of shape substrate.subdomain_nb_grid_pts or

substrate.topography_subdomain_nb_grid_pts

(default None)

tolfloat

(default None) tolerance for termination. For detailed control, use solver-specific options.

callbackcallable

(default None) callback function to be at each iteration

as callback(disp_k) where disp_k is the current displacement vector. Instead of a callable, it can be set to ‘True’, in which case the system’s default callback function is called.

logger :

(default None) log information at every objective evaluation.

property nb_grid_pts: For systems, nb_grid_pts can become non-trivial

abstract objective(offset, disp0=None, gradient=False, logger=None)

This helper method exposes a scipy.optimize-friendly interface to the evaluate() method. Use this for optimization purposes, it makes sure that the shape of disp is maintained and lets you set the offset and ‘forces’ flag without using scipy’s cumbersome argument passing interface. Returns a function of only disp Keyword Arguments: offset – determines indentation depth disp0 – preexisting displacement. influences e.g., the

physical_sizes of the proxy system in some ‘smart’ system subclasses

gradient – (default False) whether the gradient is supposed to be: used

logger – (default None) log information at every iteration.

shape_minimisation_input(in_array): For minimisation of smart systems, the initial guess array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues). Use the output of this function as argument x0 for scipy minimisation functions. Also, if you initial guess has a shape that makes no sense, this will tell you before you get caught in debugging scipy-code

shape_minimisation_output(in_array)

For minimisation of smart systems, the output array (e.g. displacement) may have a non-intuitive shape and physical_sizes (The problem physical_sizes may be decreased, as for free, non-periodic systems, or increased as with augmented-lagrangian-type issues). Use this function to get the array shape you expect to have

Arguments: in_array – array with the initial guess. has the intuitive shape you

think it has

ContactMechanics.test_surface_topography module

Module contents

Defines all solid mechanics model used in ContactMechanics