Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_integral_operators.py
2c694d3bf38ffeatures/netcdf-mpi

Last Change
Raw File

test_integral_operators.py
View Options

# -*- coding: utf-8 -*-
# @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 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
@pytest.fixture(scope="module",
params=[tm.integration_method.linear,
tm.integration_method.cutoff])
def integral_fixture(tamaas_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)
engine = model.getIntegralOperator('kelvin_force')
set_integration_method(engine, 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.getDisplacement()
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)
engine.apply(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
h = 0.01 * E # Hardening modulus
sigma_y = 0.3 * E # Initial yield stress
# Creation of model
model = tm.ModelFactory.createModel(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
residual = tm.ModelFactory.createResidual(model, sigma_y, h)
residual.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.to_voigt(eigenstress.reshape(discretization + [9]))
# Array references
gradient = residual.getVector()
stress = residual.getStress()
# Applying operator
# mindlin = model.getIntegralOperator("mindlin")
mindlin_gradient = model.getIntegralOperator("mindlin_gradient")
# Not testing displacements yet
# mindlin.apply(eigenstress, model.getDisplacement())
# Applying gradient
mindlin_gradient.apply(eigenstress, gradient)
model.applyElasticity(stress, gradient)
stress -= eigenstress
# Normalizing stess as in Mura (1976)
T, alpha = 1., 1.
beta = alpha * T * (1+nu) / (1-nu)
stress *= 1. / (2 * mu * beta)
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