Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/tests/test_integral_operators.py
787782e1adeafeatures/declarative

Last Change
Raw File

test_integral_operators.py
View Options

#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @author Lucas Frérot <lucas.frerot@epfl.ch>
#
# @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 <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
import tamaas as tm
import numpy as np
from numpy.linalg import norm
from register_integral_operators import *
def test_kelvin_volume_force():
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])]
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
c4science · Help