Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/examples/legacy/hertz.py
74b9233dc365features/netcdf-mpi

Last Change
Raw File

hertz.py
View Options

#!/usr/bin/env python
# @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/>.
import tamaas
def dumpHertzPrediction(file, bem, pressure, r, E ):
A = tamaas.SurfaceStatistics.computeContactArea(bem.getTractions())
Aratio = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions())
radius = np.sqrt(A/np.pi)
#L = bem.getSurface().getL()
L = 1.
load = pressure * L*L
radius_hertz = pow(0.75*load*r/E,1./3.)
p0_hertz = 1./np.pi*pow(6.*E*E*load/r/r,1./3.)
p0 = tamaas.SurfaceStatistics.computeMaximum(bem.getTractions())
n = bem.getTractions().shape[0]
computed_load = bem.getTractions().sum()/n/n*L*L
file.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\t{9}\n".format(pressure,load,computed_load,Aratio,A,radius,radius_hertz,radius_hertz/radius,p0,p0_hertz))
################################################################
def main(argv):
parser = argparse.ArgumentParser(description='Hertz test for the Boundary element method of Stanley and Kato')
parser.add_argument("--N",type=int, help="Surface size", required=True)
parser.add_argument("--r",type=float, help="radius of hertz sphere", required=True)
parser.add_argument("--steps",type=int, help="number of steps within which the pressure is applied.", required=True)
parser.add_argument("--precision",type=float, help="relative precision, convergence if | (f_i - f_{i-1})/f_i | < Precision.", required=True)
parser.add_argument("--e_star",type=float, help="effective elastic modulus" , required=True)
parser.add_argument("--L",type=float, help="size of the surface" , required=True)
parser.add_argument("--max_pressure",type=float, help="maximal load requested" , required=True)
parser.add_argument("--plot_surface",type=bool, help="request output of text files containing the contact pressures on the surface" , default=False)
parser.add_argument("--nthreads",type=int, help="request a number of threads to use via openmp to compute" , default=1)
args = parser.parse_args()
arguments = vars(args)
n = arguments["N"]
r = arguments["r"]
max_pressure = arguments["max_pressure"]
Ninc = arguments["steps"]
epsilon = arguments["precision"]
Estar = arguments["e_star"]
L = arguments["L"]
plot_surface = arguments["plot_surface"]
nthreads = arguments["nthreads"]
pressure = 0.
dp = max_pressure/float(Ninc)
s = np.zeros((n,n))
print s.shape
for i in range(0,n):
for j in range(0,n):
x = 1.*i - n/2
y = 1.*j - n/2
d = (x*x + y*y)*1./n/n*L*L
if d < r*r: s[i,j] = - r + np.sqrt(r*r-d)
else: s[i,j] = - r
print "\n::: DATA ::::::::::::::::::::::::::::::\n"
print " [N] {0}\n".format(n)
print " [r] {0}\n".format(r)
print " [Pext] {0}\n".format(max_pressure)
print " [Steps] {0}\n".format(Ninc)
print " [Precision] {0}\n".format(epsilon)
bem = tamaas.BemPolonski(s)
bem.setEffectiveModulus(Estar)
file = open("hertz-prediction",'w')
file.write("#pressure\tload\tcomputed-load\tarea_ratio\tarea\tradius\tradius-hertz\tradius/radius-hertz\tp0\tp0-hertz\n")
for i in range(0,Ninc):
pressure += dp
bem.computeEquilibrium(epsilon,pressure)
A = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions())
dumpHertzPrediction(file,bem,pressure,r,Estar)
if A == 1.0:
file.close()
break
tamaas.dumpTimes()
################################################################
import sys
import argparse
import numpy as np
main(sys.argv)
c4science · Help