diff --git a/language_bindings/python/muSpectre/tools.py b/language_bindings/python/muSpectre/tools.py index 9bf49ab..e2fed79 100644 --- a/language_bindings/python/muSpectre/tools.py +++ b/language_bindings/python/muSpectre/tools.py @@ -1,372 +1,375 @@ #!/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 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. + Method to compute the real space deformed structure 'final_positions' from + the deformation gradient 'F'. 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] + F=[i,j,x,y,z]; F.shape=[3,3,Nx,Ny,Nz] + where Nx, Ny, Nz are the number of gridpoints in each direction. For 1D + and 2D the form is analog. resolutions: list Grid resolutions in the Cartesian directions. For a three dimensional problem you need the following input form: - [nx, ny, nz] + [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) + TODO + ---- --> 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 - - # seperate F in homogeneous F_homo and non-homogeneous (locally fluctuating - # displacement gradient) F_nonhomo=Grad(w(x)) 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 - ### Initialize the undeformed positions "X" and the q-vectors "q" - X = np.zeros((ndim, Nx, Ny, Nz)) - q = np.zeros((ndim, Nx, Ny, Nz)) + F = np.array(F) #deformation gradient + res = np.array(resolutions) #grid resolution array + lens = np.array(lengths) #box lengths array + ndim = len(res) #dimension of the problem + d = lens/res #gridspacing in each direction - 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) + ### Initialize the undeformed positions "X" and the q-vectors "q" + # Compute a well dimensioned X and q vector (available for 1D-3D) + q_1d = lambda i: 2*np.pi*np.fft.fftfreq(res[i], d[i]) + if ndim == 1: + X = np.mgrid[0:lens[0]:d[0]] + qx = q_1d(0) + q = np.array(np.meshgrid(qx, indexing='ij')) + elif ndim == 2: + X = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1]] + qx = q_1d(0) + qy = q_1d(1) + q = np.array(np.meshgrid(qx,qy, indexing='ij')) + elif ndim == 3: + X = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1], 0:lens[2]:d[2]] + qx = q_1d(0) + qy = q_1d(1) + qz = q_1d(2) + q = np.array(np.meshgrid(qx, qy, qz, indexing='ij')) + + # Compute absolute norm square of q, |q|^2, and correct |q^2| for q=0. + # Necessary because we will compute q/|q^2|, to prevent problems we set + # |q^2|='inf' for q=0. Later the term for q=0 will be corrected. + q_norm_square = (q**2).sum(axis=0) + q_norm_square.itemset(0,float("inf")) #correct for q=0 - 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]]) - # 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") + ### Handle the deformation gradient + # Separate F in homogeneous F_homo (average deformation gradient ) and + # non-homogeneous (locally fluctuating displacement gradient) + # F_nonhomo=Grad(w(x)) deformation gradient. + 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 = np.einsum('ij,k...->ijk...', F_homo, np.ones(res)) + F_nonhomo = F-F_homo - # compute Fourier transformed of the deformation gradient, F_q = DFT(F) - F_q = np.fft.fftn(F_nonhomo, [Nx,Ny,Nz]) - F_av = F_homo + # Compute Fourier transformed of the non-homogeneous deformation gradient + # along the grid axes (s=res); F_q = DFT(F_nonhomo) + F_q = np.fft.fftn(F_nonhomo, s=res) - ### Corrections to the Fourier transformed of the derivative ### + ### 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. - # + # Method took from: # A. Vidyasagar et al., Computer Methods in Applied Mechanics and # Engineering 106 (2017) 133-151, sec. 3.4 and Appendix C. + # + # mixed with: + # P. Eisenlohr et al., International Journal of Plasticity 46 (2013) + # 37-53, Appendix B, eq. 23. 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 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)) + #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) ))) + #u_q = (F_q / (1j * (8*np.sin(q*dx)/(6*dx) - np.sin(2*q*dx)/(6*dx) ))) else: + #default order = 0 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') + 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])) - - # calculate U = *X + u - fx = np.einsum("ij...,j...->i...", F_av, X) #fx = *X - U = fx - u + nonaffine_displacements = np.real(np.fft.ifftn(u_q, s=res)) - #give more meaningfull names to the results - nonaffine_displacements = u - final_positions = U + # calculate final_positions = F_homo*X + nonaffine_displacements, + # homogeneous deformation fx = F_homo*X + fx = np.einsum("ij...,j...->i...", F_homo, X) + final_positions = fx - nonaffine_displacements 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. + computing each edge of a volume (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 50b248e..51a2566 100644 --- a/tests/python_muSpectre_tools_test.py +++ b/tests/python_muSpectre_tools_test.py @@ -1,203 +1,230 @@ #!/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 µ 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: - ----------- + Parameters + ---------- res: list, int - List of grid resoultion in each direction [Nx, Ny, Nz] + List of grid resoultion in each direction [Nx, Ny, Nz] lens: list, float - List of box lengths in each direction [Lx, Ly, Lz] + List of box lengths in each direction [Lx, Ly, Lz] - Returns: - -------- + Returns + ------- ndim: int - Dimension of the structure, derived from len(res). + Dimension of the structure, derived from len(res). X: np.ndarray, float - initial grid point positions. + initial grid point positions. F: np.ndarray - initial deformation gradient (right shape, all values set to zero) + initial deformation gradient (right shape, all values set to zero) Chi_X: np.ndarray - initial deformation field (right shape, all values set to zero) + 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 - d = lens/res #gridspacing in each direction + res = np.array(res) #grid resolution + lens = np.array(lens) #box lengths + ndim = len(res) #dimension of the problem + d = lens/res #gridspacing in each direction #setup grid point positions/initial positions X if ndim == 1: - X = np.mgrid[0:lens[0]:d[0]] + X = np.mgrid[0:lens[0]:d[0]].reshape(1,res[0]) elif ndim == 2: X = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1]] elif ndim == 3: X = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1], 0:lens[2]:d[2]] F = np.zeros((ndim,) + X.shape) Chi_X = np.zeros(X.shape) return ndim, X, F, Chi_X class MuSpectre_tools_Check(unittest.TestCase): """ Check the implementation of all muSpectre.tools functions. """ def setUp(self): self.norm_tol = 1e-8 #tolerance for the norm(Chi_X - integral(F)) def test_compute_displacements_integration(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 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 and non-diagonal deformation gradients. - Todo: - ----- + 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 """ ### 1 diagonal deformation gradient res = [8, 15, 13] lens = [4, 2.3, 2.1] 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] nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - f_pos), self.norm_tol) ### cosine diagonal deformation gradient res = [4, 2, 2] lens = [4, 2, 2] ndim, X, F, Chi_X = init_X_F_Chi(res, lens) for i in range(ndim): F[i,i,:,:,:] = 2 +0.2*2*np.pi/lens[i]*np.cos(2*np.pi*X[i]/lens[i]) Chi_X[i,:,:,:] = 2*X[i] + 0.2 * np.sin(2*np.pi*X[i]/lens[i]) nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - f_pos), 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 = [2000, 1, 1] lens = [4, 4, 4] 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,:,:,:] = 16*X[i]**2 -8*X[i]**3 +X[i]**4 if i == 0: #You have to subtract the average displacement field. This is #necessary because the Chi_X = F_av*X + w(X) and the boundary #conditions require w^-(X) = w^+(X) Chi_X[i,:,:,:] -= 128/15 nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - f_pos), self.norm_tol) ### non-diagonal deformation gradient # This one needs a lot of grid points, think more than 1000^3, therefore # I should first improve the implementation of compute_diplacements. # Before it takes to much time/crashes #res = [100, 500, 500] #lens = [8, 8, 8] #ndim, X, F, Chi_X = init_X_F_Chi(res, lens) #F[0,0,:,:,:] = 8 -6*X[0] + 3/4*X[0]**2 #F[1,1,:,:,:] = -6*X[1] + 3/4*X[1]**2 #F[2,2,:,:,:] = -6*X[2] + 3/4*X[2]**2 #F[1,0,:,:,:] = 8 #F[2,0,:,:,:] = 8 #for i in range(ndim): # Chi_X[i,:,:,:] = 8*X[0] -3*X[i]**2 +1/4*X[i]**3 # #nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) ### non-diagonal deformation gradient + ### 1D ### + res = [5] + lens = [8] + ndim, X, F, Chi_X = init_X_F_Chi(res, lens) + print(ndim, X, F, Chi_X, X.shape, F.shape, Chi_X.shape) + F[0,0,:] = 5 + 4*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]) + + nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) + + ### 2D ### + res = [5, 5] + lens = [8, 8] + ndim, X, F, Chi_X = init_X_F_Chi(res, lens) + print(ndim, X, F, Chi_X, X.shape, F.shape, Chi_X.shape) + 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]) + + nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) + + ### 3D ### res = [5, 5, 5] lens = [8, 8, 8] 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]) nonaff_dis, f_pos = µ.tools.compute_displacements(F, res, lens, order=0) #print('\nX:\n', X) #print('\nF:\n', F) #print('\nChi_X:\n', Chi_X) #print('\nfinal positions:\n', f_pos) #print('\nnonaffine displacements:\n', nonaff_dis) #print('\ndifference:\n', Chi_X - f_pos) #print('\ndifference final norm: ', np.linalg.norm(Chi_X - f_pos)) self.assertLess(np.linalg.norm(Chi_X - f_pos), self.norm_tol) if __name__ == '__main__': unittest.main()