diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py index 4241b8f..68c67d5 100644 --- a/tests/test_integral_operators.py +++ b/tests/test_integral_operators.py @@ -1,225 +1,229 @@ # -*- 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 . 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(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