Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/tests/test_integral_operators.py
6a058ccd143cmaxwell_viscoelastic…

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
import matplotlib.pyplot as plt
from pyevtk.hl import gridToVTK
def dumpGrid(name, domain, field, func=lambda x: x):
discretization = np.zeros(3, dtype=np.uint32)
discretization[0] = field.shape[0]
discretization[1] = field.shape[1] + 1
discretization[2] = field.shape[2] + 1
field_vtk = [np.zeros([discretization[i] for i in range(3)]) for j in range(3)]
coords = [np.linspace(0, domain[i], discretization[i]) for i in reversed(range(3))]
coords = np.meshgrid(*coords, indexing='ij')
for i in range(3):
field_vtk[i][:, :discretization[1]-1, :discretization[2]-1] = field[:, :, :, i]
field_vtk[i][:, discretization[1]-1, :discretization[2]-1] = field[:, 0, :, i]
field_vtk[i][:, :discretization[1]-1, discretization[2]-1] = field[:, :, 0, i]
field_vtk[i][:, discretization[1]-1, discretization[2]-1] = field[:, 0, 0, i]
field_vtk[i] = np.ascontiguousarray(func(field_vtk[i][:, :, :]))
gridToVTK(name,
*coords,
pointData={
'X':field_vtk[0],
'Y':field_vtk[1],
'Z':field_vtk[2],
})
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
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')
model = tm.ModelFactory.createModel(tm.model_type.volume_2d,
domain,
discretization)
model.E = E
model.nu = nu
displacement = model.getDisplacement()
engine = tm._tamaas._test_features.Kelvin(model)
source = np.zeros_like(displacement)
source[N//2, :, :, 2] = np.sin(omega[0]*x) * np.sin(omega[1]*y)
# source[N//2, N//2, N//2, 2] = 1.
print("Starting computation")
engine.applyVolumeForcePotential(source, displacement)
print("End computation")
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)
dumpGrid('kelvin', domain, displacement, lambda x: np.einsum('ijk->jki', x))
dumpGrid('source', domain, source, lambda x: np.einsum('ijk->jki', x))
dumpGrid('solution', domain, solution, lambda x: np.einsum('ijk->jki', x))
c4science · Help