diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py index d839764..7783219 100644 --- a/tests/test_integral_operators.py +++ b/tests/test_integral_operators.py @@ -1,72 +1,206 @@ -#!/usr/bin/env python -# coding: utf-8 -# ----------------------------------------------------------------------------- -# @author Lucas Frérot -# -# @section LICENSE -# -# Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de -# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des -# Solides) -# -# Tamaas is free software: you can redistribute it and/or modify it under the -# terms of the GNU Lesser General Public License as published by the Free -# Software Foundation, either version 3 of the License, or (at your option) any -# later version. -# -# Tamaas 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 Lesser General Public License for more -# details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with Tamaas. If not, see . -# ----------------------------------------------------------------------------- +""" +@author Lucas Frérot + +@section LICENSE + +Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de +Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des +Solides) + +Tamaas is free software: you can redistribute it and/or modify it under the +terms of the GNU Lesser General Public License as published by the Free +Software Foundation, either version 3 of the License, or (at your option) any +later version. + +Tamaas 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 Lesser General Public License for more +details. + +You should have received a copy of the GNU Lesser General Public License +along with Tamaas. If not, see . +""" + import tamaas as tm import numpy as np + from numpy.linalg import norm -from register_integral_operators import * +from register_integral_operators import register_kelvin_force + -def test_kelvin_volume_force(): +def test_kelvin_volume_force(tamaas_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') - coords = [np.linspace(0, domain[i], discretization[i], endpoint=False) for i in range(2)]\ - +[np.linspace(0, domain[2], discretization[2])] + coords = [np.linspace(0, domain[i], + discretization[i], endpoint=False) 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) + 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(tamaas_fixture): + # Definition of modeled domain + model_type = tm.model_type.volume_2d + discretization = [126, 128, 128] + system_size = [1., 3., 3.] + + # 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) + + # Coordinates + x = np.linspace(0, system_size[1], discretization[1], endpoint=False) + y = np.linspace(0, system_size[2], discretization[2], endpoint=False) + z = np.linspace(0, system_size[0], discretization[0], endpoint=True) + 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.reshape(residual.getPlasticStrain(), + discretization + [3, 3]) + eigenstress[ball, ...] = beta + eigenstress = 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]) + + 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 + ) + + z_error = norm(vertical_stress[:, 8] - 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