Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/tests/test_integral_operators.py
e28cf92a57d0maxwell_viscoelastic…

Last Change
Raw File

test_integral_operators.py
View Options

# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
# Copyright (©) 2020-2023 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/>.
import tamaas as tm
import numpy as np
import pytest
from numpy.linalg import norm
from register_integral_operators import register_kelvin_force, \
set_integration_method
def test_python_register_operator():
class Id(tm.IntegralOperator):
"""Custom Python sub-class of IntegralOperator"""
def __init__(self, model):
super().__init__(model)
self.sum = 0
def apply(self, input, output):
self.sum = int(input.sum())
output[0] = input
m = tm.Model(tm.model_type.volume_2d, [1, 1, 1], [8, 8, 8])
m.traction[:] = 1
# Creating subclass
operator = Id(m)
# Applying operator
operator(m.traction, m.displacement)
np.testing.assert_allclose(m.traction, m.displacement[0], rtol=1e-12)
assert operator.sum == 192
# Registering operator
m.operators["test"] = operator
m.traction[..., 2] += .5
# Accessing operator
m.operators["test"](m.traction, m.displacement)
assert m.operators["test"].sum == 224
@pytest.fixture(params=[tm.integration_method.linear,
tm.integration_method.cutoff],
ids=['linear', 'cutoff'])
def integral_fixture(request):
return request.param
def test_kelvin_volume_force(integral_fixture):
N = 65
E = 1.
nu = 0.3
mu = E / (2*(1+nu))
domain = np.array([1.] * 3)
omega = 2 * np.pi * np.array([1, 1]) / domain[:2]
omega_ = norm(omega)
discretization = [N] * 3
model = tm.ModelFactory.createModel(tm.model_type.volume_2d,
domain,
discretization)
model.E = E
model.nu = nu
register_kelvin_force(model)
kelvin = model.operators['kelvin_force']
set_integration_method(kelvin, integral_fixture, 1e-12)
coords = [np.linspace(0, domain[i],
discretization[i], endpoint=False, dtype=tm.dtype)
for i in range(2)]\
+ [np.linspace(0, domain[2], discretization[2])]
x, y = np.meshgrid(*coords[:2], indexing='ij')
displacement = model['displacement']
source = np.zeros_like(displacement)
# The integral of forces should stay constant
source[N//2, :, :, 2] = np.sin(omega[0]*x) * np.sin(omega[1]*y) * (N-1)
kelvin(source, displacement)
z = coords[2] - 0.5
z, x, y = np.meshgrid(z, *coords[:2], indexing='ij')
solution = np.zeros_like(source)
solution[:, :, :, 0] = -np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \
* omega[0]*z*np.cos(omega[0]*x)*np.sin(omega[1]*y)
solution[:, :, :, 1] = -np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \
* omega[1]*z*np.sin(omega[0]*x)*np.cos(omega[1]*y)
solution[:, :, :, 2] = np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \
* (3-4*nu + omega_*np.abs(z))*np.sin(omega[0]*x)*np.sin(omega[1]*y)
error = norm(displacement - solution) / norm(solution)
assert error < 5e-2
def test_mindlin(integral_fixture):
# Definition of modeled domain
# tm.set_log_level(tm.LogLevel.debug)
model_type = tm.model_type.volume_2d
discretization = [126, 128, 128]
system_size = [1., 3., 3.]
integration_method = integral_fixture
print(integration_method)
# Material contants
E = 1. # Young's modulus
nu = 0.3 # Poisson's ratio
# Creation of model
model = tm.Model(model_type, system_size, discretization)
model.E = E
model.nu = nu
mu = E / (2*(1+nu))
lamda = E * nu / ((1+nu) * (1-2*nu))
# Setup for integral operators
tm.ModelFactory.registerVolumeOperators(model)
model.setIntegrationMethod(integration_method, 1e-12)
# Coordinates
x = np.linspace(0, system_size[1], discretization[1],
endpoint=False, dtype=tm.dtype)
y = np.linspace(0, system_size[2], discretization[2],
endpoint=False, dtype=tm.dtype)
z = np.linspace(0, system_size[0], discretization[0],
endpoint=True, dtype=tm.dtype)
z, x, y = np.meshgrid(z, x, y, indexing='ij')
# Inclusion definition
a, c = 0.1, 0.2
center = [system_size[1] / 2, system_size[2] / 2, c]
r = np.sqrt((x-center[0])**2 + (y-center[1])**2 + (z-center[2])**2)
ball = r < a
# Eigenstrain definition
alpha = 1 # constant isotropic strain
beta = (3 * lamda + 2 * mu) * alpha * np.eye(3)
eigenstress = np.zeros(discretization + [3, 3], dtype=tm.dtype)
eigenstress[ball, ...] = beta
eigenstress = tm.compute.to_voigt(eigenstress.reshape(discretization + [9]))
# Array references
stress = np.zeros(discretization + [6], dtype=tm.dtype)
gradient = np.zeros_like(stress)
# Applying operator
# mindlin = model.operators["mindlin"]
mindlin_gradient = model.operators["mindlin_gradient"]
# Not testing displacements yet
# mindlin(eigenstress, model['displacement'])
# Applying gradient
mindlin_gradient(eigenstress, gradient)
model.operators['hooke'](gradient, stress)
stress -= eigenstress
# Normalizing stess as in Mura (1976)
T, alpha = 1., 1.
beta = alpha * T * (1+nu) / (1-nu)
stress *= 1. / (2 * mu * beta)
# Testing free surface
for comp in (2, 3, 4):
np.testing.assert_allclose(stress[0, ..., comp], 0, atol=1e-15)
n = discretization[1] // 2
vertical_stress = stress[:, n, n, :]
z_all = np.linspace(0, 1, discretization[0], dtype=tm.dtype)
sigma_z = np.zeros_like(z_all)
sigma_t = np.zeros_like(z_all)
inclusion = np.abs(z_all - c) < a
# Computing stresses for exterior points
z = z_all[~inclusion]
R_1 = np.abs(z - c)
R_2 = np.abs(z + c)
sigma_z[~inclusion] = 2*mu*beta*a**3/3 * (
1 / R_1**3
- 1 / R_2**3
- 18*z*(z+c) / R_2**5
+ 3*(z+c)**2 / R_2**5
- 3*(z-c)**2 / R_1**5
+ 30*z*(z+c)**3 / R_2**7
)
sigma_t[~inclusion] = 2*mu*beta*a**3/3 * (
1 / R_1**3
+ (3-8*nu) / R_2**3
- 6*z*(z+c) / R_2**5
+ 12*nu*(z+c)**2 / R_2**5
)
# Computing stresses for interior points
z = z_all[inclusion]
R_1 = np.abs(z - c)
R_2 = np.abs(z + c)
sigma_z[inclusion] = 2*mu*beta*a**3/3 * (
- 2 / a**3
- 1 / R_2**3
- 18*z*(z+c) / R_2**5
+ 3*(z+c)**2 / R_2**5
+ 30*z*(z+c)**3 / R_2**7
)
sigma_t[inclusion] = 2*mu*beta*a**3/3 * (
- 2 / a**3
+ (3-8*nu) / R_2**3
- 6*z*(z+c) / R_2**5
+ 12*nu*(z+c)**2 / R_2**5
)
# This test can be used to debug
if False:
import matplotlib.pyplot as plt
plt.plot(z_all, vertical_stress[:, 2])
plt.plot(z_all, sigma_z / (2 * mu * beta))
plt.figure()
plt.plot(z_all, vertical_stress[:, 0])
plt.plot(z_all, sigma_t / (2 * mu * beta))
plt.show()
z_error = norm(vertical_stress[:, 2] - sigma_z / (2 * mu * beta)) \
/ discretization[0]
t_error = norm(vertical_stress[:, 0] - sigma_t / (2 * mu * beta)) \
/ discretization[0]
assert z_error < 1e-3 and t_error < 1e-3
c4science · Help