diff --git a/language_bindings/python/muSpectre/tools.py b/language_bindings/python/muSpectre/tools.py index ea73824..14d95ba 100644 --- a/language_bindings/python/muSpectre/tools.py +++ b/language_bindings/python/muSpectre/tools.py @@ -1,391 +1,428 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file tools.py @author Richard Leute @date 20 Apr 2018 @brief Toolbox for computing and plotting results of muSpectre computations. @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 numpy as np import itertools +def own_fft(a, axes, x, q): + """ + Discrete Fourier transformation where it is possible to controll the x and + k-vectors. Necessary because there could be a problem... + returns out_ijk = Sum a_ijk * exp(-1j * x_ijk * q_ijk) + + Parameters: + ----------- + a: np.ndarray + vector which should be fourier transformed. + s: sequence of ints + along the last len(s) axes a will be transformed. + x: np.ndarray + x-vector which is used for the fourier transformation + q: np.ndarray + q-vector which is used for the fourier transformation. + + Return: + ------- + out: np.ndarray complex + forier transformed of a along given axes s. + """ + + #for i in s: + + out = a + + return 'under construction' + def compute_displacements(F, resolutions, lengths, order=0): """ Method to compute the real space deformed structure from the deformation gradient 'F' and the lengths 'lengths' of the system. Parameters ---------- F: list Deformation gradient. For a three dimensional problem you need the following input form: F=[i,j,x,y,z]; F.shape=[3,3,nx,ny,nz] resolutions: list Grid resolutions in the Cartesian directions. For a three dimensional problem you need the following input form: [nx, ny, nz] lengths: list Physical size of the cell in the Cartesian directions. For a three dimensional problem you need the following input form: [Lx, Ly, Lz] order: integer Order of the corrections to the Fourier transformed of the derivative. This is necessary because of finite grid spacing and hence errors in the computation of the finite difference derivative. Corrections are done as in: A. Vidyasagar et al., Computer Methods in Applied Mechanics and Engineering 106 (2017) 133-151, sec. 3.4 and Appendix C. available: order = {0, 1, 2} Returns ------- nonaffine_displacements: list Nonaffine displacement vectors 'u' pointing from the affine deformed initial positions 'X + *X' to the final position 'U' of each gridpoint. For a uniform material the nonaffine displacements should become zero. To get the nonaffine displacements 'u' one has to integrate F=du/dX by a fourier transformation. final_positions: list Positions 'U' after applying the deformation, according to the deformation gradient F, to the initial positions 'X'. U=u+*X. Caution! -------- Up to now only for 3D grids implemented! TODO: ----- --> implementation for 2D grid, (1D grid) --> implement higher order corrections, is this necessary? --> we should implement that one can give the initial coordinates 'X', if the body is already deformed. --> implementation of a test!!! """ #for 3D grids... F = np.array(F) #deformation gradient res = np.array(resolutions) lens = np.array(lengths) ndim = len(res) #dimension of the problem Nx, Ny, Nz = res #number of gridpoints in each direction Lx, Ly, Lz = lens #lengths in each direction dx, dy, dz = lens/res #gridspacing in each direction - #only part of the deformation gradient which is not one - #delta_F = F - np.eye(ndim) + print('dx, dy, dz:', dx, dy, dz) - print('\ndx/y/z: ', dx, dy, dz) + #seperate F in homogeneous F_homo and + #non-homogeneous F_nonhomo deformation gradient + F_homo = 1./(np.prod(res)) * F.sum(axis=(-1,-2,-3)) + F_homo = np.einsum('ij,k...->ijk...', F_homo, np.ones(res)) + F_nonhomo = F-F_homo + print('\nF_homo:\n', F_homo) + print('\nF_nonhomo:\n', F_nonhomo) ### Initialize the undeformed positions "X" and the q-vectors "q" X = np.zeros((ndim, Nx, Ny, Nz)) q = np.zeros((ndim, Nx, Ny, Nz)) - #own construction - qx_1d = (2*np.pi)*np.arange(Nx)/(dx*Nx) - qy_1d = (2*np.pi)*np.arange(Ny)/(dy*Ny) - qz_1d = (2*np.pi)*np.arange(Nz)/(dz*Nz) - - print('own qx_1d: ', qx_1d) - qx_1d = (np.fft.fftfreq(Nx, d=dx))*(2*np.pi) qy_1d = (np.fft.fftfreq(Ny, d=dy))*(2*np.pi) qz_1d = (np.fft.fftfreq(Nz, d=dz))*(2*np.pi) - print('np.fft qx_1d: ', qx_1d) + #qx_1d = (np.fft.fftfreq(Nx, d=1))#*(2*np.pi) + #qy_1d = (np.fft.fftfreq(Ny, d=1))#*(2*np.pi) + #qz_1d = (np.fft.fftfreq(Nz, d=1))#*(2*np.pi) + + #qx_1d = np.arange(Nx)/(dx*Nx)#*(2*np.pi) + #qy_1d = np.arange(Ny)/(dy*Ny)#*(2*np.pi) + #qz_1d = np.arange(Nz)/(dz*Nz)#*(2*np.pi) + + print('qx_1d: ', qx_1d) + print('qy_1d: ', qy_1d) + print('qz_1d: ', qz_1d) for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)): X[:,i,j,k] = np.array([i*dx, j*dy, k*dz]) q[:,i,j,k] = np.array([qx_1d[i], qy_1d[j], qz_1d[k]]) + #print('\nq:\n', q) + # compute absolute norm square of q, |q|^2, and correct for q=0 q_norm_square = (q**2).sum(axis=0) q_norm_square[0,0,0] = float("inf") # compute Fourier transformed of the deformation gradient, F_q = DFT(F) - F_q = np.fft.fftn(F, [Nx,Ny,Nz]) #should I normalize this some how?????? - F_av = np.real(1.0/(Nx*Ny*Nz) * F_q[:,:,0,0,0]) + # norm='ortho in FFT makes normalization by 1/sqrt(N)' + F_q = np.fft.fftn(F_nonhomo, [Nx,Ny,Nz]) + print('\nF_q:\n', F_q) + F_av = F_homo#np.real((Nx*Ny*Nz)**(-1.0) * F_q[:,:,0,0,0]) #1.0/(Nx*Ny*Nz) * - #delta_F_q = np.fft.fftn(delta_F, [Nx,Ny,Nz]) - #delta_F_av = np.real(1.0/(Nx*Ny*Nz) * delta_F_q[:,:,0,0,0]) - - print('average deformation gradient:\n', F_av) + #F_q = np.fft.fftn(delta_F, [Nx,Ny,Nz]) ### Corrections to the Fourier transformed of the derivative ### # Necessary because of finite grid spacing and hence errors in the # computation of the finite difference derivative. # # A. Vidyasagar et al., Computer Methods in Applied Mechanics and # Engineering 106 (2017) 133-151, sec. 3.4 and Appendix C. with np.errstate(divide='ignore', invalid='ignore'): if order == 0: # zeroth order: # computed as in: # P. Eisenlohr et al., International Journal of Plasticity 46 (2013) # 37-53, Appendix B, eq. 23. - u_q = np.einsum("ij...,j...->i...",F_q, 1j*q)/q_norm_square - #delta_u_q = np.einsum("ij...,j...->i...",delta_F_q, 1j*q)/q_norm_square + #print('\niq/q^2:\n', 1j*q/q_norm_square) + u_q = np.einsum("ij...,j...->i...", F_q, 1j*q)/q_norm_square elif order == 1: # first oder correction # correct this bei Eisenlohr print('\nThis order has to be corrected, please use order=0\n') u_q = (F_q / (1j * np.sin(q*dx)/dx)) elif order == 2: # second oder correction # correct this bei Eisenlohr print('\nThis order has to be corrected, please use order=0\n') u_q = (F_q / (1j * (8*np.sin(q*dx)/(6*dx) - np.sin(2*q*dx)/(6*dx) ))) else: print('\n\n WARNING! \n') print('In muSpectre/language_bindings/python/muSpectre/tools.py') print('The order {} is not supported. Up to now only order={{0,1,2}} is' ' supported.'.format(order)) print('I fall back to the default order: order = 0\n\n') u_q = np.einsum("ij...,j...->i...",F_q, 1j*q)/q_norm_square - # get u by inverse fourier transform u_q; u=IDFT(u_q) - u = np.real(np.fft.ifftn(u_q, [Nx,Ny,Nz])) - #delta_u = np.real(np.fft.ifftn(delta_u_q, [Nx,Ny,Nz])) + print('\nu_q:\n', u_q, u_q.shape) - print('\nu:\n', u) + # get u by inverse fourier transform u_q; u=IDFT(u_q) + u = np.real(np.fft.ifftn(u_q, [Nx,Ny,Nz])) # (2*np.pi)**(-ndim/2) * # calculate U = *X + u fx = np.einsum("ij...,j...->i...", F_av, X) #fx = *X #delta_fx = np.einsum("ij...,j...->i...", delta_F_av, X) #I think this is nonsense print('\nfx:\n', fx) #print('\ndelta_fx:\n', delta_fx) U = fx - u #give more meaningfull names to the results nonaffine_displacements = u final_positions = U return nonaffine_displacements, final_positions def point_to_volume_grid(point_grid, resolutions, lengths, affine_displacements=np.array(None)): """ Function which converts the point grid (grid on which computations are done) into a volume grid (for visualization). The interpolation is done by computing each edge of a volume, which center is the gridpoint, from the average of the eight (for 3D) surrounding gridpoints. For the boundaries the periodicity of the problem is used. Parameters ---------- point_grid: Real np.ndarray (with dimension of the problem) grid around which the volume cells are constructed. The problem is assumed to be periodic. resolutions: list Grid resolutions in the Cartesian directions. For a three dimensional problem you need the following input form: [nx, ny, nz] lengths: list Physical size of the cell in the Cartesian directions. For a three dimensional problem you need the following input form: [Lx, Ly, Lz] affine_displacements: Real np.ndarray (with dimension of the problem) gives the affine displacements on each point of the point_grid. They are used to fix the average over affine displaced points at the boundaries. Returns ------- volume_grid: Real np.ndarray (with dimension of the problem) nodes of all cells constructed around the points of the point_grid. Therefore the array has in each dimension one entry more than the point_grid. TODO ---- --> improve this for 2D grids and write a general code --> shorten the loop expressions """ if affine_displacements.all() == None: #fix them, by setting all to zero (no affine deformation of the grid) affine_displacements = np.zeros(point_grid.shape) # weighted sum of point_grid; # The position of an edge is given by: 1/#NN * Sum_{i element NN edge} r_i . # Thus the weighted (1/#NN) positions r_i of all nearest points of the # point_grid. dim = len(resolutions) Lx, Ly, Lz = lengths Nx, Ny, Nz = resolutions summ = np.zeros(point_grid.shape) weight = 0 ### compute the edges of each volume as the average position of the 8(3D) or # 4(2D) surrounding 'point_grid' points. If the point grid originates from # an affine + a non-affine deformation the affine deformation has to be # subtracted/added at the boundaries. for d in range(dim+1): for ax in itertools.combinations(range(1,dim+1),d): #Move one of the neighbouring 8 points on the central point rolled_p_g = np.roll(point_grid, 1, axis=ax) #Corresponding affine deformation to the rolled point_grid rolled_aff_d = np.roll(affine_displacements, 1, axis=ax) #Difference between the rolled and unrolled affine deformations. #This is used for the correction at the periodic boundaries. aff_d_difference = affine_displacements - rolled_aff_d #change sign and value of the first entry of the rolled axis because #this is the pbc entry and thus is shifted by the lengths in the #direction of the axis. #Substract/add also the affine part of the deformation. You need #only the affine deformation planar to the pbc surface, because #perpendicular the displacement is given by the box size L. for i in ax: if i == 1: rolled_p_g[0 , 0,:,:] -= Lx #rolled_p_g[: , 0,:,:] += aff_d_difference[:, 0,:,:] #rolled_p_g[1: , 0,:,:] -= rolled_aff_d[1:, 0,:,:] #rolled_p_g[1: , 0,:,:] -= affine_displacements[1:, 0,:,:] elif i == 2: rolled_p_g[1 , :,0,:] -= Ly #rolled_p_g[: , :,0,:] += aff_d_difference[:, :,0,:] #rolled_p_g[0 , :,0,:] -= rolled_aff_d[0, :,0,:] #rolled_p_g[2 , :,0,:] -= rolled_aff_d[2, :,0,:] #rolled_p_g[0 , :,0,:] -= affine_displacements[0, :,0,:] #rolled_p_g[2 , :,0,:] -= affine_displacements[2, :,0,:] elif i == 3: rolled_p_g[2 , :,:,0] -= Lz #rolled_p_g[: , :,:,0] += aff_d_difference[:, :,:,0] #rolled_p_g[:2 , :,:,0] -= rolled_aff_d[:2, :,:,0] #rolled_p_g[:2 , :,:,0] -= affine_displacements[:2, :,:,0] summ += rolled_p_g weight += 1 volume_grid_no_pbc = 1/weight * summ ### extend the volume_grid by the periodic images of the boundary points #shape of the periodic grid(extended in each direction by 1) new_shape = np.array(point_grid.shape) new_shape[1:] += 1 new_shape = new_shape.tolist() volume_grid = np.zeros(new_shape) #fill the grid by the old values volume_grid[:, :Nx, :Ny, :Nz] = volume_grid_no_pbc #Fill the grid by the periodic images. #And extend the affine displacements by one point in each direction. #The first displacement is the periodic one to the the last displacement in #each direction. affine_displacements_ext = np.zeros(new_shape) affine_displacements_ext[:, :Nx, :Ny, :Nz] = affine_displacements for d in range(dim): if d == 0: values = np.copy(volume_grid[:, 0, :Ny, :Nz]) affine_displacements_ext[:, -1, :Ny, :Nz] = \ affine_displacements[:, 0, :Ny, :Nz] values[0] += Lx #values[1:] += affine_displacements_ext[1:, 0, :Ny, :Nz] volume_grid[:, -1, :Ny, :Nz] = values elif d == 1: values = np.copy(volume_grid[:, :, 0, :Nz]) affine_displacements_ext[:, :Nx, -1, :Nz] = \ affine_displacements[:, :, 0, :Nz] values[1] += Ly #values[0] += affine_displacements_ext[0, :, 0, :Nz] #values[2] += affine_displacements_ext[2, :, 0, :Nz] volume_grid[:, :, -1, :Nz] = values elif d == 2: values = np.copy(volume_grid[:, :, :, 0]) affine_displacements_ext[:, :Nx, :Ny, -1] = \ affine_displacements[:,:,:,0] values[2] += 0#Lz #values[:2] += affine_displacements_ext[:2,:,:,0] volume_grid[:, :, :, -1] = values return volume_grid def write_structure_as_vtk(file_name, positions, cellData=None, pointData=None): """ Function to save a structure as vtk files by using pyevtk.hl. Two files a written named: file_name_grid.vts and file_name_points.vtu. Parameters ---------- file_name: str name of the 'vtk' file which will be writen positions: list grid point positions cellData: dictionary data associated with voxel properties pointData: dictionary data associated to each point/point properties Returns ------- writes a '.vts' file with the file name, file_name_grid.vts and a '.vtu' file with the name file_name_points.vtu. If it ends successfully it returns the stirng 'file_name_grid.vts and file_name_points.vtu files were written.' TODO: --> improve this... --> implement 2D and 3D """ from pyevtk.hl import gridToVTK, pointsToVTK dim = len(positions.shape)-1 print('dimension: ', dim) cdata = cellData pdata = pointData if dim >= 2: x = positions[0] y = positions[1] # set all z positions to zero z = None #z = positions[0] #z[:] = 0.0 if dim == 3: z = positions[2] #print(x.shape, y.shape, z.shape) if dim == 2: gridToVTK(file_name+'_grid', x, y, cellData=cdata, pointData=pdata) pointsToVTK(file_name+'_points', x, y, data=pdata) elif dim == 3: gridToVTK(file_name+'_grid', x, y, z, cellData=cdata, pointData=pdata) pointsToVTK(file_name+'_points', x, y, z, data=pdata) message = '{}_grid.vts and {}_points.vtu'\ ' files were written.'.format(file_name, file_name) return message diff --git a/tests/python_muSpectre_tools_test.py b/tests/python_muSpectre_tools_test.py index c2f48e0..5018e58 100644 --- a/tests/python_muSpectre_tools_test.py +++ b/tests/python_muSpectre_tools_test.py @@ -1,125 +1,132 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file python_muSpectre_tools_test.py @author Richard Leute @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 itertools from python_test_imports import µ class MuSpectre_tools_Check(unittest.TestCase): """ Check the implementation of all tools. """ def setUp(self): - #self.resolutions = [4, 3, 1]#[4, 3, 1],[4,5,1] - #self.lengths = [2*np.pi, 4*np.pi, 4*np.pi] - self.resolutions = [4, 1, 1] - self.lengths = [4, 2, 2] + #set for Phi_X[i,:,:,:] = 2*X[i] + 0.2 * np.sin(X[i]*(2/3)) + #self.resolutions = [4, 1, 1] + #self.lengths = [3*np.pi, 3*np.pi, 3*np.pi] + + #set for Phi_X[i,:,:,:] = X[i] + 0.5*(X[i]**2.0) - 1./3*(X[i]**3) + #self.resolutions = [5, 1, 1] + #self.lengths = [1.5, 1.5, 1.5] + + #set for Phi_X[i,:,:,:] = 16*X[i]**2 -8*X[i]**3 +X[i]**4 + self.resolutions = [2000, 1, 1] + self.lengths = [4, 4, 4] def test_compute_displacements_integration(self): """ - Proof if compute_displacements turns a cosine into a sine and a 1+x^2 - into x+x^3 for different grid resolutions and and cell lengths. - Additionally it checks the implementation of the correction order. + Proof if compute_displacements turns a deforamtion gradient F into its + corresponding deformation field Chi_X. This is tested for different grid + resolutions and 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. This is tested for 2D and 3D grids. - test off diagonal deformation gradients + test off diagonal and non-diagonal deformation gradients. + + Todo: + ----- + * Implementation of other deformation gradients which are periodic. + Why it doesn't work for e.g. 0.5*x^2-1/3*x^3 for L=1.5. At 1.5 and 0 + the deformation is zero and thus the signal can be periodically """ - # create a cosine deformation gradient + # initialization of parameters res = np.array(self.resolutions) lens = np.array(self.lengths) ndim = len(res) #dimension of the problem Nx, Ny, Nz = res #number of gridpoints in each direction Lx, Ly, Lz = lens #lengths in each direction dx, dy, dz = lens/res #gridspacing in each direction X = np.zeros((ndim, Nx, Ny, Nz)) for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)): X[:,i,j,k] = np.array([i*dx, j*dy, k*dz]) - print('\nX:\n', X) - - # compute deformation gradient F and final positions Phi_X - # F_ij = dx_j/dX_i ; x=final position, X=initial position - # x_i = Phi(X_i); Let Phi(X_i)=sin(X_i) - # ==> F_ij = dx_j/dX_i = dsin(X_j)/dX_i = cos(X_j)*delta_ij + # compute deformation fields Chi_X and corresponding deformation + # gradients F. (We stick here close to the notation as in the paper + # from P. Eisenlohr). F = np.zeros((ndim, ndim, Nx, Ny, Nz)) - F_2 = np.zeros((ndim, ndim, Nx, Ny, Nz)) - Phi_X = np.zeros((ndim, Nx, Ny, Nz)) - Phi_X_2 = np.zeros((ndim, Nx, Ny, Nz)) + Chi_X = np.zeros((ndim, Nx, Ny, Nz)) for i in range(ndim): - F[i,i,:,:,:] = np.cos(X[i]) - Phi_X[i,:,:,:] = np.sin(X[i]) + #F[i,i,:,:,:] = np.cos(X[i]) + #Phi_X[i,:,:,:] = np.sin(X[i]) - F_2[i,i,:,:,:] = 1+X[i]# + np.cos(X[i]) - Phi_X_2[i,:,:,:] = X[i]+0.5*X[i]**2# + np.sin(X[i]) + #F[i,i,:,:,:] = 1 + 0.2*np.cos(X[i]) + #Phi_X[i,:,:,:] = X[i] + 0.2*np.sin(X[i]) - #print('\nF:\n', F) - # correct Phi_X by the affine displacements - F_av = 1/(Nx*Ny*Nz)*F.sum(axis=(2,3,4)) - #Phi_final_X = np.einsum('ij...,j...->i...', F_av, X) + Phi_X #think this is nonsense - #print('\nF_av:\n', F_av) - nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) + #F[i,i,:,:,:] = 2 + 0.2*(2/3)*np.cos(X[i]*(2/3)) + #Phi_X[i,:,:,:] = 2*X[i] + 0.2 * np.sin(X[i]*(2/3)) + + #F[i,i,:,:,:] = 1 + X[i] + #Phi_X[i,:,:,:] = X[i] + 0.5*(X[i]**2.0) - F_2_av = 1/(Nx*Ny*Nz)*F_2.sum(axis=(2,3,4)) - #Phi_final_X_2 = np.einsum('ij...,j...->i...', F_2_av, X) + Phi_X_2 #think this is nonsense - nonaff_dis_2, f_pos_2 = µ.tools.compute_displacements(F_2, res, lens, order=0) + #F[i,i,:,:,:] = 1 + X[i] - X[i]**2 + #Phi_X[i,:,:,:] = X[i] + 0.5*(X[i]**2.0) - 1./3*(X[i]**3) - # compare with sin of np.linspace - test_res_x = np.sin(np.linspace(0.0, Lx, Nx, endpoint=False)) - test_res_y = np.sin(np.linspace(0.0, Ly, Ny, endpoint=False)) - test_res_z = np.sin(np.linspace(0.0, Lz, Nz, endpoint=False)) - #print('\nSin(X[0]):\n', test_res_x) - #print('\nSin(X[1]):\n', test_res_y) - #print('\nSin(X[2]):\n', test_res_z) + #F[i,i,:,:,:] = X[i] - X[i]**2 + #Phi_X[i,:,:,:] = 0.5*(X[i]**2.0) - 1./3*(X[i]**3) + + #F[i,i,:,:,:] = 1 + #Phi_X[i,:,:,:] = X[i] + + F[i,i,:,:,:] = 32*X[i] -24*X[i]**2 +4*X[i]**3 + Chi_X[i,:,:,:] = 16*X[i]**2 -8*X[i]**3 +X[i]**4 + if i == 0: + Chi_X[i,:,:,:] -= 128/15 #subtract the average displacement field + + nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) #print('\nfinal positions:\n', f_pos) #print('\nnonaffine displacements:\n', nonaff_dis) #print('\nPhi(X):\n', Phi_X) - #print('\nPhi_final(X):\n', Phi_final_X) #print('\ndifference:\n', Phi_X - f_pos) - #print('\ndifference final norm: ', np.linalg.norm(Phi_final_X - f_pos)) - - print('\n\n F_2:\n') - print('\nF_2_av:\n', F_2_av) - print('\nfinal positions 2:\n', f_pos_2) - print('\nnonaffine displacements 2:\n', nonaff_dis_2) - print('\nPhi(X) 2:\n', Phi_X_2) - #print('\nPhi_final(X) 2:\n', Phi_final_X_2) - print('\ndifference 2:\n', Phi_X_2 - f_pos_2) - print('\ndifference 2 norm: ', np.linalg.norm(Phi_X_2 - f_pos_2)) - #print('\ndifference 2 final norm: ', np.linalg.norm(Phi_final_X_2 - f_pos_2)) - - self.assertLess(np.linalg.norm(Phi_X - f_pos), 1e-8) + #print('\ndifference final norm: ', np.linalg.norm(Phi_X - f_pos)) + + self.assertLess(np.linalg.norm(Chi_X - f_pos), 1e-8) if __name__ == '__main__': unittest.main()