Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/tests/test_matvec.py
5722d08fe801maxwell_viscoelastic…

Last Change
Raw File

test_matvec.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 pytest
import numpy as np
import tamaas as tm
from numpy.testing import assert_allclose
from test_dumper import model_fixture # noqa
linalg = pytest.importorskip('scipy.sparse.linalg')
def test_matvec_westergaard(model_fixture):
m = model_fixture
m.be_engine.registerNeumann()
op = m.operators['Westergaard::neumann']
linop = linalg.aslinearoperator(op)
x = [np.linspace(0, L, N, endpoint=False)
for (L, N) in zip(m.boundary_system_size, m.boundary_shape)]
coords = np.meshgrid(*x, indexing='ij')
t = m.traction
res = np.zeros_like(m.displacement)
volume = len(m.shape) != len(m.boundary_shape)
basic = m.type in (tm.model_type.basic_1d, tm.model_type.basic_2d)
# Adapting to special cases
if basic:
t = t[..., np.newaxis]
res = res[..., np.newaxis]
t[:] = 0
t[..., -1] = np.sin(2 * np.pi * coords[0])
# Computing using native API
op(t, res)
# Only look at surface displacement
if volume:
res = res[0]
# Computing using LinearOperator API
y = np.reshape(linop * np.ravel(t), res.shape)
assert_allclose(res, y)
def test_matvec_hooke():
m = tm.Model(tm.model_type.volume_2d, [1, 1, 1], [2, 2, 2])
op = m.operators['hooke']
linop = linalg.aslinearoperator(op)
strain = np.zeros(m.shape + [6], dtype=tm.dtype)
strain[..., 3:] = 1
stress = np.zeros_like(strain)
op(strain, stress)
res = linop * np.ravel(strain)
assert_allclose(res, np.ravel(stress))
assert_allclose(res, np.ravel(strain) * m.mu * 2)
def test_matvec_mindlin():
m = tm.Model(tm.model_type.volume_2d, [1, 1, 1], [2, 2, 2])
tm.ModelFactory.registerVolumeOperators(m)
op = m.operators['mindlin_gradient']
linop = linalg.aslinearoperator(op)
strain = np.zeros(m.shape + [6], dtype=tm.dtype)
strain[..., 3:] = 1
stress = np.zeros_like(strain)
op(strain, stress)
res = linop * np.ravel(strain)
assert_allclose(res, np.ravel(stress))
c4science · Help