Diffusion tamaas (features/netcdf-mpi) Browse

netcdf-mpi

adhesion.py
View Options

			#!/usr/bin/env python
			# @file
			# @section LICENSE
			#
			# Copyright (©) 2016-2020 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/>.

			import sys
			import tamaas as tm
			import numpy as np
			import matplotlib.pyplot as plt

			def plotSurface(surface):
			fig = plt.figure()
			ax = fig.add_subplot(111)
			img = ax.imshow(surface)
			#fig.colorbar(img)

			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_size2 - (x2 + y**2))
			return disp.copy()

			def main():
			grid_size = 128
			curvature = 0.5
			effective_modulus = 1.
			load = 0.1

			surface_energy = 0
			rho = 2.071e-7

			surface = constructHertzProfile(grid_size, curvature)
			# SG = tm.SurfaceGeneratorFilterFFT()
			# SG.getGridSize().assign(grid_size)
			# SG.getHurst().assign(0.8)
			# SG.getRMS().assign(0.002);
			# SG.getQ0().assign(8);
			# SG.getQ1().assign(8);
			# SG.getQ2().assign(16);
			# SG.getRandomSeed().assign(156);
			# SG.Init()
			# surface = SG.buildSurface()
			print "Max height {}".format(surface.max())
			print "Min height {}".format(surface.min())

			bem = tm.BemGigi(surface)
			bem.setDumpFreq(1)
			functional = tm.ExponentialAdhesionFunctional(bem)
			functional.setParameter('rho', rho)
			functional.setParameter('surface_energy', surface_energy)
			bem.setEffectiveModulus(effective_modulus)
			bem.addFunctional(functional)
			bem.computeEquilibrium(1e-6, load)

			tractions = bem.getTractions()
			print "Average pressure = {}".format(tractions.mean())
			# bem.computeTrueDisplacements()
			t_displacements = bem.getTrueDisplacements()
			t_gap = bem.getGap()

			plotSurface(tractions)

			plt.figure()
			plt.plot(surface[grid_size/2, :])
			plt.title("Surface")
			plt.figure()
			plt.plot(tractions[grid_size/2, :])
			plt.title("Pressure")
			plt.figure()
			plt.plot(t_gap[grid_size/2, :])
			plt.title("Gap")
			plt.figure()
			plt.plot(t_displacements[grid_size/2, :])
			plt.title("Displacement")
			plt.figure()
			plt.plot(t_displacements[grid_size/2, :] - surface[grid_size/2, :])
			plt.title("Disp-surf")
			plotSurface(t_displacements)

			plt.show()
			return 0

			# Testing contact area against Hertz solution for solids of revolution
			contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
			hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
			hertz_area = np.pi * hertz_contact_size**2
			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 * load * effective_modulus*2 curvature 2)(1. / 3.) / np.pi
			pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
			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, tractions[i, j] > 0) for i in range(grid_size) for j in range(grid_size)]
			contact_indexes = map(lambda x: x[0:2], filter(lambda x: x[2], contact_indexes))

			# 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 += 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] - 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)

			if area_error > 1e-2 or pressure_error > 1e-2 or displacement_error > 1e-4:
			return 1

			return 0

			if __name__ == "__main__":
			sys.exit(main())

tamaas/examples/legacy/adhesion.py19806151094ffeatures/netcdf-mpi

adhesion.pyView Options

tamaas/examples/legacy/adhesion.py
19806151094ffeatures/netcdf-mpi

adhesion.py
View Options