Diffusion tamaas (features/netcdf-mpi) Browse

netcdf-mpi

test_hertz_disp.py
View Options

			# -- coding: utf-8 --
			# @file
			# @section LICENSE
			#
			# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
			# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
			#
			# This program is free software: you can redistribute it and/or modify
			# it under the terms of the GNU Affero General Public License as published
			# by the Free Software Foundation, either version 3 of the License, or
			# (at your option) any later version.
			#
			# This program 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 Affero General Public License for more details.
			#
			# You should have received a copy of the GNU Affero General Public License
			# along with this program. If not, see <https://www.gnu.org/licenses/>.

			from __future__ import print_function, division

			import tamaas as tm
			import numpy as np


			def constructHertzProfile(size, curvature):
			radius = 1. / curvature
			x = np.linspace(-0.5, 0.5, size)
			y = np.linspace(-0.5, 0.5, size)
			x, y = np.meshgrid(x, y)
			surface = np.sqrt(radius2 - x2 - y**2)
			surface -= surface.min()
			return surface.copy()


			def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
			x = np.linspace(-0.5, 0.5, size)
			y = np.linspace(-0.5, 0.5, size)
			x, y = np.meshgrid(x, y)
			disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 -
			(x2 + y2))
			return disp.copy()


			def test_hertz_disp():
			tm.initialize()
			grid_size = 512
			curvature = 0.5
			effective_modulus = 1.
			load = 0.12

			surface = constructHertzProfile(grid_size, curvature)
			model = tm.ModelFactory.createModel(tm.model_type_basic_2d,
			[1., 1.],
			[grid_size, grid_size])
			model.setElasticity(1, 0)
			solver = tm.PolonskyKeerRey(model, surface, 1e-12,
			tm.PolonskyKeerRey.gap,
			tm.PolonskyKeerRey.gap)

			solver.solve(load - surface.mean())

			tractions = model.getTraction()
			true_displacements = model.getDisplacement()
			true_gap = true_displacements - surface
			pressure = np.mean(tractions)
			contact_points = np.where(1.0e-10 >= true_gap)
			contact_area = (np.size(contact_points)/2)/float(grid_size*grid_size)

			print('The contact area computed with the gap map is '+str(contact_area))

			hertz_contact_size = (3 * pressure / (4 * curvature * effective_modulus))**(1. / 3.)

			print('Exact Hertz contact radius is '+str(hertz_contact_size))
			hertz_area = np.pi * hertz_contact_size**2
			print('Exact Hertz contact area is '+str(hertz_area))
			area_error = np.abs(hertz_area - contact_area) / hertz_area
			print("Area error: {}".format(area_error))

			# Testing maximum pressure
			max_pressure = tractions.max()
			hertz_max_pressure = (6 * pressure * effective_modulus*2 curvature 2)(1. / 3.) / np.pi
			pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
			print('Exact Hertz max pressure is '+str(hertz_max_pressure))
			print('max pressure is '+str(np.max(tractions)))
			print("Max pressure error: {}".format(pressure_error))

			# Testing displacements
			hertz_displacements = computeHertzDisplacement(effective_modulus,
			hertz_contact_size,
			hertz_max_pressure,
			grid_size)

			# Selecing only the points that are in contact
			contact_indexes = [(i, j)
			for i in range(grid_size) for j in range(grid_size)
			if true_gap[i, j] < 1e-12]

			# Displacements of bem are centered around the mean of the whole surface
			# and Hertz displacements are not centered, so we need to compute mean
			# on the contact zone for both arrays
			bem_mean = 0.
			hertz_mean = 0.
			for index in contact_indexes:
			bem_mean += true_displacements[index]
			hertz_mean += hertz_displacements[index]

			bem_mean /= len(contact_indexes)
			hertz_mean /= len(contact_indexes)
			# Correction applied when computing error
			correction = hertz_mean - bem_mean

			# Computation of error of displacement in contact zone
			error = 0.
			hertz_norm = 0.
			for index in contact_indexes:
			error += (hertz_displacements[index] - true_displacements[index] - correction)**2
			hertz_norm += (hertz_displacements[index] - hertz_mean)**2

			displacement_error = np.sqrt(error / hertz_norm)
			print("Displacement error (in contact zone): {}".format(displacement_error))

			assert area_error < 1e-2 and pressure_error < 1e-2 and displacement_error < 2e-3

			tm.finalize()


			if __name__ == "__main__":
			test_hertz_disp()

tamaas/tests/test_hertz_disp.pydea0fa68046bfeatures/netcdf-mpi

test_hertz_disp.pyView Options

tamaas/tests/test_hertz_disp.py
dea0fa68046bfeatures/netcdf-mpi

test_hertz_disp.py
View Options