Page MenuHomec4science
Files F60505145

test_hertz.py
No OneTemporary
Actions

File Metadata

Created
Tue, Apr 30, 17:22

test_hertz.py
View Options

# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
# Copyright (©) 2020-2024 Lucas Frérot
#
# 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 division, print_function
import numpy as np
from numpy.linalg import norm
import tamaas as tm
def compute_perimeter(tractions):
contact = np.zeros_like(tractions, bool)
contact[tractions > 0] = True
# Computing clusters
clusters = tm.FloodFill.getClusters(contact, False)
return np.sum([c.perimeter for c in clusters])
def test_contact_area(pkr):
"Testing contact area against Hertz solution for solids of revolution"
model, hertz = pkr
def error(erzt, true):
return np.abs(erzt - true) / true
tractions = model['traction']
true_area = hertz.a**2 * np.pi
contact_area = tm.Statistics2D.contact(tractions)
hertz_area = np.mean(hertz.pressure > 0)
assert error(contact_area, hertz_area) < 1e-2
assert error(contact_area, true_area) < 5e-3
contact_area = tm.Statistics2D.contact(tractions,
compute_perimeter(tractions))
hertz_perimeter = compute_perimeter(hertz.pressure)
hertz_area -= (np.pi - 1 + np.log(2)) / 24 \
* hertz_perimeter / tractions.size
assert error(contact_area, hertz_area) < 1e-2
assert error(contact_area, true_area) < 3e-3
def test_pressure(pkr):
model, hertz = pkr
error = norm(model['traction'] - hertz.pressure) / norm(hertz.pressure)
assert error < 5e-3
def test_displacement(pkr):
model, hertz = pkr
error = norm(model['displacement'] - hertz.displacement)\
/ norm(hertz.displacement)
assert error < 8e-3

Event Timeline

c4science · Help