Homec4science
Diffusion µSpectre (master) Browse

muSpectre/tests/python_muSpectre_tools_test.py
03da37d380d8master

Last Change
Raw File

python_muSpectre_tools_test.py
View Options

#!/usr/bin/env python3
# -*- coding:utf-8 -*-
"""
@file python_muSpectre_tools_test.py
@author Richard Leute <richard.leute@imtek.uni-freiburg.de>
@date 24 May 2018
@brief test the behaviour of tools.py
@section LICENSE
Copyright © 2018 Till Junge
µSpectre is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 3, or (at
your option) any later version.
µSpectre 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
General Public License for more details.
You should have received a copy of the GNU General Public License
along with GNU Emacs; see the file COPYING. If not, write to the
Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.
"""
import unittest
import numpy as np
import scipy.misc as sm
import tempfile
import itertools
import os
from python_test_imports import µ
def init_X_F_Chi(res, lens):
"""
Setup all the needed parameters for initialization of the deformation
gradient F and the corresponding deformation map/field Chi_X.
Parameters
----------
res: list, int
List of grid resoultion in each direction [Nx, Ny, Nz]
lens: list, float
List of box lengths in each direction [Lx, Ly, Lz]
Returns
-------
ndim: int
Dimension of the structure, derived from len(res).
X: np.ndarray, float
initial grid point positions.
F: np.ndarray
initial deformation gradient (right shape, all values set to zero)
Chi_X: np.ndarray
initial deformation field (right shape, all values set to zero)
"""
# initialization of parameters
res = np.array(res) #grid resolution
lens = np.array(lens) #box lengths
ndim = len(res) #dimension of the problem
X = µ.tools.compute_real_space_lattice(res, lens)
F = np.zeros((ndim,) + X.shape)
Chi_X = np.zeros(X.shape)
return ndim, X, F, Chi_X
def nth_order_central_diff_derivative(function_values, grid_spacing, order):
"""
Compute the first derivative of a function with values 'function_values'
with the central difference approximation to the order 'order'. The function
values are on a rectangualr grid with constant grid spacing 'grid_spacing'
for each direction.
CAUTION:
The function is assumed to be periodic (pbc)!
Thus, if there is a discontinuity at the boundaries you have to expect
errors in the derivative at the vicinity close to the dicontinuity.
Parameters
----------
function_values: numpy.ndarray, float
value of the function on an equally spaced grid,
with grid spacing 'grid_spacing'.
grid_spacing: float or np.array of floats
gives the grid spacing distance in each grid direction. For a single
value the grid spacing is assumed similar in each direction.
order: int, order >= 1
gives the accuracy order of the central difference approximation. This
means for first order (order=1), three grid points are used to compute
the derivative; order=2, 5 grid points; etc.
Returns
-------
deriv: np.ndarray, float
central difference derivative
"""
f = function_values
d = grid_spacing
n = order
dim = function_values.shape[0] #dimension of the problem (1D-3D)
# Get the stencils/ central difference weights
Np = 2*n + 1 #number of gridpoints used
weights = sm.central_diff_weights(Np)
# Central difference approx.
deriv = np.zeros(shape=(dim,)+f.shape)
for ax in range(dim): #get derivatives in each direction (x,y,z)
for i in range(Np): #get the fuction values on the left/right hand side
deriv[:,ax] += weights[i] * np.roll(f, n-i, axis=ax+1) * 1/d[ax]
return deriv
class MuSpectre_tools_Check(unittest.TestCase):
"""
Check the implementation of all muSpectre.tools functions.
TODO
----
* better tests for:
-> test_point_to_volume_grid() (off diagonal F and nonaff displacements)
-> test_write_structure_as_vtk()
Check not only if the files were written but also their content!
"""
def setUp(self):
self.lengths = np.array([2.4, 3.7, 4.1])
self.resolutions = np.array([8, 9, 11])
self.norm_tol = 1e-8 #tol for norm(Chi_X - integral(F))
self.norm_tol_ho = 5e-2 #tol for problematic higher order integr.
self.norm_tol_cent_diff = 1e-8 #tol for central-diff. approx. test
self.norm_tol_i_DFT = 1e-8 #tol for i_DFT test
def test_nth_order_central_diff_derivative(self):
"""
Test of the nth order central difference approximation of the first
derivative of a function on a grid.
"""
### 1D
res = np.array([50])
lens = np.array([7])
d = lens/res
ndim, x, deriv, f = init_X_F_Chi(res, lens)
f = np.sin(2*np.pi/lens * x) #f(x)=axsin(a*x)
approx_deriv = nth_order_central_diff_derivative(f, grid_spacing=d,
order=4)
deriv = 2*np.pi/lens * np.cos(2*np.pi/lens * x) #f'(x)=a*cos(a*x)
self.assertLess(np.linalg.norm(deriv-approx_deriv),
self.norm_tol_cent_diff)
### 3D
res = np.array([46, 60, 41])
lens = np.array([4, 1.7, 5.3])
d = lens/res
ndim, x, deriv, f = init_X_F_Chi(res, lens)
for i in range(ndim):
deriv[i,i,:,:,:] = 2*np.pi/lens[i] * np.cos(2*np.pi*x[i]/lens[i])
f[i,:,:,:] = np.sin(2*np.pi*x[i]/lens[i])
approx_deriv = nth_order_central_diff_derivative(f, grid_spacing=d,
order=5)
self.assertLess(np.linalg.norm(deriv-approx_deriv),
self.norm_tol_cent_diff)
def test_compute_real_space_lattice(self):
"""
Test the function compute_real_space_lattice which creates an orthogonal
equally spaced grid of the given resolution and lengths.
"""
res = [3, 30, 5]
lens = [5, 2.2, 7]
shift = - np.array([0.2, 0.3, 0.4]).reshape((3,1,1,1))
d_offset = np.array([0.1, 0.2, 0.3])
Nx, Ny, Nz = res
dx, dy, dz = np.array(lens)/np.array(res)
new_d = np.array(lens-d_offset)/(np.array(res)-1)
new_dx, new_dy, new_dz = new_d
grid = np.zeros((3,)+tuple(res)) #control grid
grid_d = np.zeros((3,)+tuple(res)) #other grid spacing
for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)):
grid[:,i,j,k] = np.array([i*dx, j*dy, k*dz])
grid_d[:,i,j,k] = np.array([i*new_dx, j*new_dy, k*new_dz])
grid_s = grid + shift
for n in range(1,4):
#loop from 1D to 3D
comp_g = µ.tools.compute_real_space_lattice(res[:n], lens[:n])
comp_g_s= µ.tools.compute_real_space_lattice(res[:n], lens[:n],
shift=shift)
comp_g_d= µ.tools.compute_real_space_lattice(res[:n], lens[:n],
d=d_offset[:n])
s = (np.s_[:n],) + (np.s_[:],)*n + (0,)*(3-n)#slice to control grids
self.assertLess(np.linalg.norm(comp_g - grid[s]), self.norm_tol)
self.assertLess(np.linalg.norm(comp_g_s - grid_s[s]), self.norm_tol)
self.assertLess(np.linalg.norm(comp_g_d - grid_d[s]), self.norm_tol)
def test_compute_reciprocal_space_lattice(self):
"""
Test the function compute_reciprocal_space_lattice which creates an
orthogonal equally spaced reciprocal grid.
"""
res = [3, 4, 5]
lens = [4, 6, 7]
Nx, Ny, Nz = res
q_1d = lambda i: (2*np.pi)/lens[i] * \
np.append(np.arange(res[i]-res[i]//2),
-np.arange(1, res[i]//2 + 1)[::-1])
qx = q_1d(0)
qy = q_1d(1)
qz = q_1d(2)
q = np.zeros((3,)+tuple(res))
for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)):
q[:,i,j,k] = np.array([qx[i], qy[j], qz[k]])
for n in range(1,4):
#loop from 1D to 3D
comp_q = µ.tools.compute_reciprocal_space_lattice(res[:n],
lens[:n])
comp_q_f = µ.tools.compute_reciprocal_space_lattice(res[:n],
lens[:n],
factor=1/2)
s = (np.s_[:n],) + (np.s_[:],)*n + (0,)*(3-n)#slice to control grids
self.assertLess(np.linalg.norm(comp_q - q[s]), self.norm_tol)
self.assertLess(np.linalg.norm(comp_q_f - q[s]/(4*np.pi)),
self.norm_tol)
def test_extend_grid_by_pbc(self):
"""
Test extend_grid_by_pbc with the function compute_real_space_lattice().
"""
res = np.array([3, 4, 5])
lens = np.array([4, 6, 7])
d = lens/res
for n in range(1,4):
grid_pbc = µ.tools.compute_real_space_lattice(res[:n], lens[:n])
#shift is lens
grid_sh = µ.tools.compute_real_space_lattice(res[:n]+1,
lens[:n]+d[:n])
grid = np.copy(grid_sh) #further modify boundaries
if n == 1:
grid[0,-1] = grid[0,0]
elif n == 2:
grid[0,-1,:] = grid[0,0,:]
grid[1,:,-1] = grid[1,:,0]
elif n == 3:
grid[0,-1,:,:] = grid[0,0,:,:]
grid[1,:,-1,:] = grid[1,:,0,:]
grid[2,:,:,-1] = grid[2,:,:,0]
comp_g_shift = µ.tools.extend_grid_by_pbc(grid_pbc, pbc_shift=lens)
comp_g = µ.tools.extend_grid_by_pbc(grid_pbc)
self.assertLess(np.linalg.norm(comp_g - grid),
self.norm_tol)
self.assertLess(np.linalg.norm(comp_g_shift - grid_sh),
self.norm_tol)
def test_i_DFT(self):
"""
Test the inverse discrete Fourier transformation (i_DFT) by comparing it
with numpy.fft.ifftn() and the given input function.
"""
for d in range(1,4):
res = self.resolutions[:d]
lens = self.lengths[:d]
ndim = len(res)
X = µ.tools.compute_real_space_lattice(res, lens)
### real signal
f_real = np.sin(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X) + 4
f_real_q = np.fft.fftn(f_real, s=res)
f_real_numpy = np.fft.ifftn(f_real_q, s=res)
f_real_reconstr = µ.tools.i_DFT(f_real_q, res, lens, X)
self.assertLess(np.linalg.norm(f_real - f_real_reconstr),
self.norm_tol_i_DFT)
self.assertLess(np.linalg.norm(f_real_numpy - f_real_reconstr),
self.norm_tol_i_DFT)
#shifted output
sh = 1/3 * lens/res #shift by 1/3 of the grid spacing
X_shift = µ.tools.compute_real_space_lattice(res, lens, shift = sh)
f_real_shift = np.sin(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) *
X_shift) + 4
f_real_shift_reconstr = µ.tools.i_DFT(f_real_q, res, lens, X_shift)
self.assertLess(np.linalg.norm(f_real_shift - f_real_shift_reconstr),
self.norm_tol_i_DFT)
### complex signal
comp = 1j * np.cos(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X)
f_comp = np.sin(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X)+ comp+2
f_comp_q = np.fft.fftn(f_comp, s=res)
f_comp_numpy = np.fft.ifftn(f_comp_q, s=res)
f_comp_reconstr = µ.tools.i_DFT(f_comp_q, res, lens, X)
self.assertLess(np.linalg.norm(f_comp - f_comp_reconstr),
self.norm_tol_i_DFT)
self.assertLess(np.linalg.norm(f_comp_numpy - f_comp_reconstr),
self.norm_tol_i_DFT)
#shifted output
comp_shift = 1j * np.cos(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) *
X_shift)
f_comp_shift = np.sin(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) *
X_shift) + comp_shift + 2
f_comp_shift_reconstr = µ.tools.i_DFT(f_comp_q, res, lens, X_shift)
self.assertLess(np.linalg.norm(f_comp_shift - f_comp_shift_reconstr),
self.norm_tol_i_DFT)
def test_compute_displacements(self):
"""
Proof if compute_displacements turns a deforamtion gradient F into its
corresponding deformation field Chi_X. This is tested for different grid
resolutions and cell lengths. One has to keep in minde that the
displacement fluctuation field w(x) has to fulfill periodicity conditions
of the form: w^-=w^+ on the surfaces delB^- , delB^+ of the periodic
region B (see: P. Eisenlohr et al., International Journal of Plasticity
46 (2013), chap. 2.1 and appendix B; The mentioned w here corresponds to
w tilde (displacement fluctuation field) in the paper).
Additionally the implementation of the correction order is checked.
Tests are done for 1D, 2D and 3D grids, with diagonal and non-diagonal
deformation gradients.
One additional tests with a highly non continuous functions is done
to check the Gibbs ringing effect. Therefore we choose a Heaviside step
function deformation gradient to integrate.
"""
order_o = [1, 2] #tested higher orders (zeroes order is always tested)
### F=1I, diagonal deformation gradient
res = [9, 15, 13]
lens = [4, 2.3, 2.1]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
for i in range(ndim):
F[i,i,:,:,:] = 1
Chi_X[i,:,:,:] = X[i]
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1,2,3)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = \
µ.tools.compute_displacements(F, res, lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
### cosine diagonal deformation gradient
### 1D ###
res = [51]
lens = [8]
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
d = np.array(lens)/np.array(res)
F[0,0,:] = 4 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
Chi_X[0,:] = 4*X + 2*np.sin(2*np.pi*X[0]/lens[0])
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = \
µ.tools.compute_displacements(F, res, lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
### non-diagonal deformation gradient
### 2D ###
res = [21, 21] #Here's a problem for odd grids because then the higer
#order solution decouples and oscillates(see plot below)
lens = [7, 1.4]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
for i in range(ndim):
F[i,i,:,:] = 2+ 0.2*4*np.pi/lens[i] * np.cos(4*np.pi*X[i]/lens[i])
Chi_X[i,:,:] = 2*X[i] + 0.2 * np.sin(4*np.pi*X[i]/lens[i])
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1,2)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = \
µ.tools.compute_displacements(F, res, lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
### 3D ###
res = [21, 21, 21] #Here's a problem for odd grids because the higher
#order solution decouples and oscillates
lens = [7, 1.4, 3]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
for i in range(ndim):
F[i,i,:,:,:] = 2 + 0.2*4*np.pi/lens[i] * np.cos(4*np.pi*X[i]
/lens[i])
Chi_X[i,:,:,:] = 2*X[i] + 0.2 * np.sin(4*np.pi*X[i]/lens[i])
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1,2,3)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = \
µ.tools.compute_displacements(F, res, lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
### polynomial diagonal deformation gradient
# Choose the prefactors of the polynomial such that at least Chi_X and F
# have respectively the same values at the boundaries (here X_i=0 and
# X_i=4).
res = [51, 51, 51]
lens = [4, 4, 4]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
for i in range(ndim):
F[i,i,:,:,:] = 32*X[i] -24*X[i]**2 +4*X[i]**3
Chi_X[i,:,:,:] = -128/15 + 16*X[i]**2 -8*X[i]**3 +X[i]**4
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol_ho)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
initial_pos, displacements = µ.tools.compute_displacements(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos+displacements)),
self.norm_tol_ho*0.016)
### non-diagonal deformation gradient
### 2D ###
res = [5, 5]
lens = [8, 8]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
F[0,0,:,:] = 5 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
F[1,1,:,:] = 5 + 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1])
F[1,0,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
for i in range(ndim):
Chi_X[i,:,:] = 5*X[i] + np.sin(2*np.pi*X[i]/lens[i]) \
+ np.sin(2*np.pi*X[0]/lens[0])
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1,2)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = µ.tools.compute_displacements(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
### 3D ###
res = [5, 5, 5]
lens = [8, 8, 8]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
F[0,0,:,:,:] = 5 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
F[1,1,:,:,:] = 5 + 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1])
F[2,2,:,:,:] = 5 + 2*np.pi/lens[2]*np.cos(2*np.pi/lens[2]*X[2])
F[1,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
F[2,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
for i in range(ndim):
Chi_X[i,:,:,:] = 5*X[i] + np.sin(2*np.pi*X[i]/lens[i]) \
+ np.sin(2*np.pi*X[0]/lens[0])
# zeroth order correction
initial_pos, displacements = µ.tools.compute_displacements(F, res, lens,
order=0)
self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)),
self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1,2,3)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = µ.tools.compute_displacements(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
#### REAL TEST, shear of material with two different Young moduli. ####
#tolerance for fourier integration with two materials
norm_tol_2mat_o0 = 2e-4
#tolerance for higher order integration with two materials
norm_tol_2mat_ho = 5e-3
### initialize
#structure
res = [5, 121]
lens = [10, 10]
ndim = len(res) #dimension of the problem
Nx, Ny = res
Lx, Ly = lens
gs = np.array(lens) / np.array(res) #grid spacing
formulation = µ.Formulation.small_strain
Young = [10, 20] #Youngs modulus for each phase (soft,hard)
Poisson = [0.3, 0.3] #Poissons ratio for each phase
order = [0] + order_o #correction order of Fourier integration
#solver
newton_tol = 1e-6 #tolerance for newton algo
cg_tol = 1e-6 #tolerance for cg algo
equil_tol = 1e-6 #tolerance for equilibrium
maxiter = 100
verbose = 0
### geometry (two slabs above each other in y-direction)
phase = np.zeros((Nx,Ny), dtype=int) #0 for soft material
boundary = res[1]//2
phase[:, boundary:] = 1 #1 for hard material
phase = phase.flatten()
cell = µ.Cell(res, lens, formulation)
mat = µ.material.MaterialLinearElastic4_2d.make(cell, "material")
for i, pixel in enumerate(cell):
#add Young and Poisson depending on the material index
m_i = phase[i] #m_i = material index (phase index)
mat.add_pixel(pixel, Young[m_i], Poisson[m_i])
cell.initialise()
DelF = np.array([[0 , 0.01],
[0.01 , 0 ]])
### compute the total deformation gradient F and integrate it
solver_newton = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose)
r_newton = µ.solvers.newton_cg(cell, DelF, solver_newton,
newton_tol, equil_tol, verbose)
grad = r_newton.grad.T.reshape(*res, 2, 2).transpose((3,2,0,1))
one = np.eye(ndim).reshape((ndim,)*2+(1,)*ndim)
F = grad + one
fin_pos = {}
for n in order:
initial_pos, displacements = \
µ.tools.compute_displacements(F, res, lens, order=n)
fin_pos[str(n)] = initial_pos + displacements #computed positions
### analytic solution of the problem
# The analytic signal doesn't know about the periodic boundary condition,
# thus the analytic box is L_ana = L-d; and the Heaviside step is
# inbetween two grid points.
#
# solution : Fourier Analytic
# Grid points : + + + # # # o + + + # # #
# Step signal : ---------\ /- ----------|
# \ / |
# \---------/ |----------
# Length : |<--------- L --------->| |<- L_ana = L-d ->|
# soft/hard : |< L_soft >|< L_hard >| | L_soft | L_hard |
# grid spacing: |gs|
l_soft_ana = gs[1]*(boundary - 1/2) #analytic height soft
l_hard_ana = gs[1]*(Ny - boundary - 1/2) #analytic height hard
l_soft = gs[1] * boundary #periodic height soft
l_hard = gs[1] * (Ny-boundary) #periodic height hard
mu = np.array(Young) / (2 * (1+np.array(Poisson))) #shear modulus
gamma = 2*DelF[0,1] #mean deformation
g_soft = (Ly*gamma) / (l_soft + l_hard * mu[0]/mu[1])
g_hard = (Ly*gamma) / (l_soft * mu[1]/mu[0] + l_hard)
###compute d (final positions)
# initial grid
X = µ.tools.compute_real_space_lattice(res, lens)
d = np.zeros(X.shape) #initialize new positions
#soft
d[0, :, :boundary] = g_soft/2 * X[1, :, :boundary]
#hard
d[0, :, boundary:] = g_hard/2 * (X[1, :, boundary:]-l_soft_ana) \
+ g_soft/2 * l_soft_ana
#y-coordinate
d[1, :, :] = gamma/2 * X[0, :, :]
#choose the right integration constant
# shift the analytic solution such that:
# -> the average nonaffine deformation is zero (integral of the
# nonaffine deformation gradient + integration const).
dim_axis = tuple(-np.arange(1,ndim+1)) #the last ndim entries
F_homo = 1./(np.prod(res)) * F.sum(axis=dim_axis)
F_homo = F_homo.reshape((ndim,)*2+(1,)*ndim)
# integration constant = - (integral / number of gridpoints)
int_const = - ((d[0] - F_homo[0,1] * X[1]).sum(axis=1))[0] / Ny
#final positions = displacements + old positions + integration const.
ana_sol = d + X + np.array([int_const, 0]).reshape((ndim,)+(1,)*ndim)
#Visualization of the test
#import matplotlib.pyplot as plt
#plt.plot(X[1,0,:], ana_sol[0,0,:], '-*', label='x analytic')
#plt.plot(X[1,0,:], fin_pos['0'][0,0,:], '--+', label='x numeric')
#plt.plot(X[1,0,:], ana_sol[1,0,:], '-*', label='y analytic')
#plt.plot(X[1,0,:], fin_pos['0'][1,0,:], '--+', label='y numeric')
#plt.plot(X[1,0,:], F[0,1,0,:], '-.o', label='F')
#plt.plot(X[1,0,:], ana_sol[0,0,:] - fin_pos['0'][0,0,:], '-+',
# label='ana-num')
#plt.legend()
#plt.show()
#check the zero order correction
norm = np.linalg.norm(fin_pos['0'] - ana_sol)
self.assertLess(norm, norm_tol_2mat_o0)
#check the higher order corrections
for n in order[1:]:
norm = np.linalg.norm(fin_pos[str(n)] - ana_sol)
self.assertLess(norm, norm_tol_2mat_ho)
def test_compute_displacements_staggered_grid(self):
"""
Proof if compute_displacements_staggered_grid turns a deforamtion
gradient F into its corresponding INTERSTITIAL deformation field
Chi_X_inter.
This is tested for different grid resolutions and cell lengths. One has
to keep in minde that the displacement fluctuation field w(x) has to
fulfill periodicity conditions of the form: w^-=w^+ on the surfaces
delB^- , delB^+ of the periodic region B (see: P. Eisenlohr et al.,
International Journal of Plasticity 46 (2013), chap. 2.1 and appendix B;
The mentioned w here corresponds to w tilde (displacement fluctuation
field) in the paper).
Additionally the implementation of the correction order is checked.
Tests are done for 1D, 2D and 3D grids, with diagonal and non-diagonal
deformation gradients.
One additional tests with a highly non continuous functions is done
to check the Gibbs ringing effect. Therefore we choose a Heaviside step
function deformation gradient to integrate.
"""
order_o = [1,2] #tested higher orders (zeroes order is always tested)
### F=1I, diagonal deformation gradient
res = [9, 15, 13]
lens = [4, 2.3, 2.1]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
for i in range(ndim):
F[i,i,:,:,:] = 1
Chi_X[i,:,:,:] = X[i]
Chi_X_inter[i,:,:,:] = X_inter[i]
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X_inter - (Chi_X_inter.sum(axis=(1,2,3))\
/np.prod(np.array(res))
+ d/2).reshape((ndim,)+(1,)*ndim)
#at the left boundary the step is integrated to zero,
#the average between left and right boundary.
Chi_X_zero_mean[0,0,:,:] = 0
Chi_X_zero_mean[1,:,0,:] = 0
Chi_X_zero_mean[2,:,:,0] = 0
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
if n == 1:
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol)
if n== 2:
#necessary because step at boundaries leads to ringing
self.assertLess(np.linalg.norm(Chi_X_zero_mean
- (initial_pos + displacements)), self.norm_tol_ho*1e2)
### cosine diagonal deformation gradient
### 1D ###
res = [27]
lens = [8]
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
d = np.array(lens)/np.array(res)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
F[0,0,:] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
Chi_X[0,:] = 2*np.sin(2*np.pi*X[0]/lens[0])
Chi_X_inter[0,:] = 2*np.sin(2*np.pi*X_inter[0]/lens[0])
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)) , self.norm_tol_ho)
### 2D ###
res = [81, 81]
lens = [7, 1.4]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
for i in range(ndim):
F[i,i,:,:] = 0.2*4*np.pi/lens[i] * np.cos(4*np.pi*X[i]/lens[i])
Chi_X[i,:,:] = 0.2 * np.sin(4*np.pi*X[i]/lens[i])
Chi_X_inter[i,:,:] = 0.2 * np.sin(4*np.pi*X_inter[i]/lens[i])
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol_ho)
### 3D ###
res = [23, 21, 28]
lens = [7, 1.4, 3]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
for i in range(ndim):
F[i,i,:,:,:] = 0.2*4*np.pi/lens[i] * np.cos(4*np.pi*X[i]/lens[i])
Chi_X[i,:,:,:] = 0.2 * np.sin(4*np.pi*X[i]/lens[i])
Chi_X_inter[i,:,:,:] = 0.2 * np.sin(4*np.pi*X_inter[i]/lens[i])
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
Chi_X_zero_mean = Chi_X - Chi_X.sum(axis=(1,2,3)).reshape((ndim,)
+(1,)*ndim)\
/np.prod(np.array(res))
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol_ho*21)
### polynomial diagonal deformation gradient
# Choose the prefactors of the polynomial such that at least Chi_X and F
# have respectively the same values at the boundaries (here X_i=0 and
# X_i=4).
res = [51, 51, 51]
lens = [4, 4, 4]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
for i in range(ndim):
F[i,i,:,:,:] = 32*X[i] - 24*X[i]**2 + 4*X[i]**3
Chi_X[i,:,:,:] = -128/15 + 16*X[i]**2 - 8*X[i]**3 + X[i]**4
Chi_X_inter[i,:,:,:] = -128/15 + 16*X_inter[i]**2 \
- 8*X_inter[i]**3 + X_inter[i]**4
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol_ho*1.6)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
#a way to large grid is needed to reduce the error to 0.05
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol_ho*140)
### non-diagonal deformation gradient
### 2D ###
res = [35, 32]
lens = [8, 8]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
F[0,0,:,:] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
F[1,1,:,:] = 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1])
F[1,0,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
for i in range(ndim):
Chi_X[i,:,:] = np.sin(2*np.pi*X[i]/lens[i]) \
+ np.sin(2*np.pi*X[0]/lens[0])
Chi_X_inter[i,:,:] = np.sin(2*np.pi*X_inter[i]/lens[i]) \
+ np.sin(2*np.pi*X_inter[0]/lens[0])
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol_ho*5)
### 3D ###
res = [41, 41, 32]
lens = [8, 8, 8]
d = np.array(lens)/np.array(res)
ndim, X, F, Chi_X = init_X_F_Chi(res, lens)
X_inter = µ.tools.compute_real_space_lattice(res, lens, shift=-d/2)
Chi_X_inter = np.zeros((ndim,)+tuple(res))
F[0,0,:,:,:] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
F[1,1,:,:,:] = 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*X[1])
F[2,2,:,:,:] = 2*np.pi/lens[2]*np.cos(2*np.pi/lens[2]*X[2])
F[1,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
F[2,0,:,:,:] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0])
for i in range(ndim):
Chi_X[i,:,:,:] = np.sin(2*np.pi*X[i]/lens[i]) \
+ np.sin(2*np.pi*X[0]/lens[0])
Chi_X_inter[i,:,:,:] = np.sin(2*np.pi*X_inter[i]/lens[i]) \
+ np.sin(2*np.pi*X_inter[0]/lens[0])
# zeroth order correction
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res, lens, order=0)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol)
# higher order correction
for n in order_o:
F = nth_order_central_diff_derivative(Chi_X, d, order=n)
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
self.assertLess(np.linalg.norm(Chi_X_inter -
(initial_pos + displacements)), self.norm_tol_ho*30)
#### REAL TEST, shear of material with two different Young moduli. ####
#tolerance for fourier integration with two materials
norm_tol_2mat_o0 = 6e-4
#tolerance for higher order integration with two materials (key=order)
norm_tol_2mat_ho = {'1': 1e-8, '2': 5e-4}
### initialize
#structure
res = [5, 21]
lens = [10, 10]
ndim = len(res) #dimension of the problem
Nx, Ny = res
Lx, Ly = lens
gs = np.array(lens) / np.array(res) #grid spacing
formulation = µ.Formulation.small_strain
Young = [10, 20] #Youngs modulus for each phase (soft,hard)
Poisson = [0.3, 0.3] #Poissons ratio for each phase
order = [0] + order_o #correction order of Fourier integration
#solver
newton_tol = 1e-6 #tolerance for newton algo
cg_tol = 1e-6 #tolerance for cg algo
equil_tol = 1e-6 #tolerance for equilibrium
maxiter = 100
verbose = 0
### geometry (two slabs above each other in y-direction)
phase = np.zeros((Nx,Ny), dtype=int) #0 for soft material
boundary = res[1]//2
phase[:, boundary:] = 1 #1 for hard material
phase = phase.flatten()
cell = µ.Cell(res, lens, formulation)
mat = µ.material.MaterialLinearElastic4_2d.make(cell, "material")
for i, pixel in enumerate(cell):
#add Young and Poisson depending on the material index
m_i = phase[i] #m_i = material index (phase index)
mat.add_pixel(pixel, Young[m_i], Poisson[m_i])
cell.initialise()
DelF = np.array([[0 , 0.01],
[0.01 , 0 ]])
### compute the total deformation gradient F and integrate it
solver_newton = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose)
r_newton = µ.solvers.newton_cg(cell, DelF, solver_newton,
newton_tol, equil_tol, verbose)
grad = r_newton.grad.T.reshape(*res, 2, 2).transpose((3,2,0,1))
one = np.eye(ndim).reshape((ndim,)*2+(1,)*ndim)
F = grad + one
fin_pos = {}
for n in order:
initial_pos, displacements = \
µ.tools.compute_displacements_staggered_grid(F, res,
lens, order=n)
fin_pos[str(n)] = initial_pos + displacements #computed positions
### analytic solution of the problem
# The analytic signal doesn't know about the periodic boundary condition,
# thus the analytic box is L_ana = L-d; and the Heaviside step is
# inbetween two grid points.
#
# solution : Fourier Analytic
# Grid points : + + + # # # o + + + # # #
# Step signal : ---------\ /- ----------|
# \ / |
# \---------/ |----------
# Length : |<--------- L --------->| |<- L_ana = L-d ->|
# soft/hard : |< L_soft >|< L_hard >| | L_soft | L_hard |
# grid spacing: |gs|
l_soft_ana = gs[1]*(boundary - 1/2) #analytic height soft
l_hard_ana = gs[1]*(Ny - boundary - 1/2) #analytic height hard
l_soft = gs[1] * boundary #periodic height soft
l_hard = gs[1] * (Ny-boundary) #periodic height hard
mu = np.array(Young) / (2 * (1+np.array(Poisson))) #shear modulus
gamma = 2*DelF[0,1] #mean deformation
g_soft = (Ly*gamma) / (l_soft + l_hard * mu[0]/mu[1])
g_hard = (Ly*gamma) / (l_soft * mu[1]/mu[0] + l_hard)
###compute d (final positions)
# initial grid
X = µ.tools.compute_real_space_lattice(res, lens, shift=-gs/2)
d = np.zeros(X.shape) #initialize new positions
#soft
d[0, :, :boundary] = g_soft/2 * X[1, :, :boundary]
#hard
d[0, :, boundary:] = g_hard/2 * (X[1, :, boundary:]-l_soft_ana) \
+ g_soft/2 * l_soft_ana
#y-coordinate
d[1, :, :] = gamma/2 * X[0, :, :]
#choose the right integration constant
# shift the analytic solution such that:
# -> the average nonaffine deformation is zero (integral of the
# nonaffine deformation gradient + integration const).
dim_axis = tuple(-np.arange(1,ndim+1)) #the last ndim entries
F_homo = 1./(np.prod(res)) * F.sum(axis=dim_axis)
F_homo = F_homo.reshape((ndim,)*2+(1,)*ndim)
# integration constant = - (integral / number of gridpoints)
int_const = - ((d[0] - F_homo[0,1] * X[1]).sum(axis=1))[0] / Ny
#final positions = displacements + old positions + integration const.
ana_sol = d + X + np.array([int_const, 0]).reshape((ndim,)+(1,)*ndim)
#Visualization of the test
#import matplotlib.pyplot as plt
#o = '2' #order of the visualized correction
#plt.plot(X[1,1,:], ana_sol[0,1,:], '-*', label='x analytic')
#plt.plot(X[1,1,:], fin_pos[o][0,1,:], '--+', label='x numeric')
#plt.plot(X[1,1,:], ana_sol[1,1,:], '-*', label='y analytic')
#plt.plot(X[1,1,:], fin_pos[o][1,1,:], '--+', label='y numeric')
#plt.plot(X[1,1,:], F[0,1,1,:], '-.o', label='F')
#plt.plot(X[1,1,:], ana_sol[0,1,:] - fin_pos[o][0,1,:], '-+',
# label='ana-num')
#plt.legend()
#plt.show()
#check the zero order correction
norm = np.linalg.norm(fin_pos['0'] - ana_sol)
self.assertLess(norm, norm_tol_2mat_o0)
#check the higher order corrections
for n in order[1:]:
norm = np.linalg.norm(fin_pos[str(n)] - ana_sol)
self.assertLess(norm, norm_tol_2mat_ho[str(n)])
def test_point_to_volume_grid(self):
"""
Test the implementation of the function point_to_volume_grid()
TODO
----
*check affine displacements/and first implement them right...
--> Do I need the nonaffine_displacement correction any more???!!!
"""
lens = self.lengths
res = self.resolutions
ndim = len(res)
d = lens/res
point_grid = np.mgrid[0:lens[0]:d[0],
0:lens[1]:d[1],
0:lens[2]:d[2]]
# for an orthogonal same spaced point_grid one can easily compute the
# volume_grid by just increase (pbc=False) the mesh by one in each
# direction and choose the points inbetween (d/2) the point_grid
# positions.
vg = np.mgrid[0-d[0]/2:lens[0]:d[0],
0-d[1]/2:lens[1]:d[1],
0-d[2]/2:lens[2]:d[2]]
vg_pbc = np.mgrid[0-d[0]/2:lens[0]-d[0]/2:d[0],
0-d[1]/2:lens[1]-d[1]/2:d[1],
0-d[2]/2:lens[2]-d[2]/2:d[2]]
# loop over 1D to 3D, for pbc=False/True
for dim in range(1,4):
s = list(np.s_[:dim,:,:,:])
try:
s[dim+1] = 0
s[dim+2] = 0
except IndexError:
pass
s = tuple(s)
F_eye = np.einsum('ij,k...->ijk...',np.eye(dim),np.ones(res[:dim]))
calc_vg = µ.tools.point_to_volume_grid(
point_grid = point_grid[s],
resolutions = res[:dim],
lengths = lens[:dim],
F = F_eye,
pbc = False,
nonaffine_displacements = np.array(None) )
self.assertLess(np.linalg.norm(calc_vg - vg[s]), self.norm_tol)
calc_vg_pbc = µ.tools.point_to_volume_grid(
point_grid = point_grid[s],
resolutions = res[:dim],
lengths = lens[:dim],
F = F_eye,
pbc = True,
nonaffine_displacements = np.array(None) )
self.assertLess(np.linalg.norm(calc_vg_pbc - vg_pbc[s]),
self.norm_tol )
def test_point_to_volume_grid_fourier(self):
"""
Test the implementation of the function point_to_volume_grid_fourier().
Should do the same as point_to_volume_grid().
"""
lens = self.lengths
res = self.resolutions
ndim = len(res)
d = lens/res
point_grid = np.mgrid[0:lens[0]:d[0],
0:lens[1]:d[1],
0:lens[2]:d[2]]
# for an orthogonal same spaced point_grid one can easily compute the
# volume_grid by just increase (pbc=False) the mesh by one in each
# direction and choose the points inbetween (d/2) the point_grid
# positions.
vg = np.mgrid[0-d[0]/2:lens[0]:d[0],
0-d[1]/2:lens[1]:d[1],
0-d[2]/2:lens[2]:d[2]]
vg_pbc = np.mgrid[0-d[0]/2:lens[0]-d[0]/2:d[0],
0-d[1]/2:lens[1]-d[1]/2:d[1],
0-d[2]/2:lens[2]-d[2]/2:d[2]]
# loop over 1D to 3D, for pbc=False/True
for dim in range(1,4):
s = list(np.s_[:dim,:,:,:])
try:
s[dim+1] = 0
s[dim+2] = 0
except IndexError:
pass
s = tuple(s)
F_eye = np.einsum('ij,k...->ijk...',np.eye(dim),np.ones(res[:dim]))
calc_vg = µ.tools.point_to_volume_grid_fourier(
point_grid = point_grid[s],
resolutions = res[:dim],
lengths = lens[:dim],
F = F_eye,
pbc = False )
self.assertLess(np.linalg.norm(calc_vg - vg[s]), self.norm_tol)
calc_vg_pbc = µ.tools.point_to_volume_grid_fourier(
point_grid = point_grid[s],
resolutions = res[:dim],
lengths = lens[:dim],
F = F_eye,
pbc = True )
self.assertLess(np.linalg.norm(calc_vg_pbc - vg_pbc[s]),
self.norm_tol )
def test_write_structure_as_vtk(self):
"""
Test the implementation of the function write_structure_as_vtk().
It will write 1D, 2D and 3D structures in a folder which will be deleted
after the files were selfuccessfully written.
TODO
----
Check not only if the files were written but also their content!
"""
#initialize some structure
lens = self.lengths
res = self.resolutions
pos = {} #positions dictionary
cd = {} #cell data dictionary
pd = {} #point data dictionary
for i in range(1,4):
name = "{}d".format(i)
pos[name] = init_X_F_Chi(res[:i], lens[:i])[1]
cd_shape = tuple(np.array(pos[name].shape[1:])-1) + (1,)*(3-i)
pd_shape = pos[name].shape[1:] + (1,)*(3-i)
cd[name] = {'random values' : np.random.random_sample(cd_shape)}
pd[name] = {'random values' : np.random.random_sample(pd_shape)}
#create a temporary directory and run the tests in it
dir_name = 'write_structure_as_vtk_tests'
cwd_start = os.getcwd()
#The temporary directory is atomatically cleand up after one is exiting
#the with block.
with tempfile.TemporaryDirectory(dir=cwd_start) as dir_name:
#move to tmp dir
os.chdir(dir_name)
#start with the test (write file and check if they were written)
for i in range(1,4):
name = "{}d".format(i)
µ.tools.write_structure_as_vtk('{}'.format(name),
pos[name], None, None)
µ.tools.write_structure_as_vtk('{}_cdpd'.format(name),
pos[name], cd[name], pd[name])
#test if files were written
error_message = '\nThere is a problem in the function ' + \
'muSpectre.tools.write_structure_as_vtk().\n' +\
'The file {} was not written.'
assert os.path.exists('{}.vts'.format(name)) == 1, \
error_message.format(name+'_grid.vts')
assert os.path.exists('{}_cdpd.vts'.format(name)) == 1, \
error_message.format(name+'_cdpd_grid.vts')
os.chdir(cwd_start)
if __name__ == '__main__':
unittest.main()
c4science · Help