Page MenuHomec4science
Files F91584235

test_bem_grid.py
No OneTemporary
Actions

File Metadata

Created
Tue, Nov 12, 11:01

test_bem_grid.py
View Options

# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# 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/>.
from __future__ import print_function
import tamaas as tm
import numpy as np
from numpy.linalg import norm
def make_profile(f, profile, k):
x = np.linspace(0, 1, profile.shape[0], endpoint=False, dtype=tm.dtype)
y = np.linspace(0, 1, profile.shape[1], endpoint=False, dtype=tm.dtype)
x, y = np.meshgrid(x, y)
profile[:, :] = f(2*k*np.pi*y)
def sin_profile(profile, k):
make_profile(np.sin, profile, k)
def cos_profile(profile, k):
make_profile(np.cos, profile, k)
def test_bem_grid(tamaas_fixture):
"""Checking that all the pipes are clean"""
size = 128
k = 1
E = 0.91
nu = 0.3
mu = E / (2*(1+nu))
pressure = np.zeros((size, size), dtype=tm.dtype)
surface = np.zeros((size, size), dtype=tm.dtype)
westergaard_disp = np.zeros((size, size), dtype=tm.dtype)
cos_profile(pressure, k)
cos_profile(surface, k)
bem = tm.BemGridPolonski(surface)
bem.setElasticity(E, nu)
bem.computeInfluence()
bem.getTractions()[:, :, 2] = pressure[:, :]
bem.computeDisplacementsFromTractions()
displacements = bem.getDisplacements()
# Constructing vertical displacement solution
cos_profile(westergaard_disp, k)
westergaard_disp *= (1-nu**2)/(E*np.pi*k)
westergaard_norm = (1-nu**2)/(E*np.pi*k)*np.sqrt(0.5) * size
# Computing error
error = norm(displacements[:, :, 2] - westergaard_disp)/westergaard_norm
assert error < 1e-10, "Error in normal displacement"
# Constructing horizontal displacement solution
sin_profile(westergaard_disp, k)
westergaard_disp *= (2*nu-1)/(4*np.pi*k*mu)
westergaard_norm = np.abs(2*nu-1)/(4*np.pi*k*mu)*np.sqrt(0.5)*size
error = norm(displacements[:, :, 1] - westergaard_disp)/westergaard_norm
assert error < 1e-10, "Error in horizontal displacement"
# Looking at normal vectors
normals = np.zeros((size, size, 3), dtype=tm.dtype)
normals_norm = np.zeros((size, size, 3), dtype=tm.dtype)
sin_profile(normals[:, :, 1], k)
normals[:, :, 1] *= 2*np.pi*k
normals[:, :, 2] = 1
x = np.linspace(0, 1, size, endpoint=False, dtype=tm.dtype)
y = np.linspace(0, 1, size, endpoint=False, dtype=tm.dtype)
x, y = np.meshgrid(x, y)
normals_norm[:, :, 0] = np.sqrt(1+4*np.pi**2*k**2*np.sin(2*np.pi*k*y)**2)
normals_norm[:, :, 1] = np.sqrt(1+4*np.pi**2*k**2*np.sin(2*np.pi*k*y)**2)
normals_norm[:, :, 2] = np.sqrt(1+4*np.pi**2*k**2*np.sin(2*np.pi*k*y)**2)
normals /= normals_norm
bem.computeSurfaceNormals()
bem_normals = bem.getSurfaceNormals()
error = norm(bem_normals - normals)
assert error < 1e-9, "Error in surface normals"
if __name__ == "__main__":
test_bem_grid()

Event Timeline

c4science · Help