diff --git a/language_bindings/python/muSpectre/tools.py b/language_bindings/python/muSpectre/tools.py index 8ff66a9..30a8d00 100644 --- a/language_bindings/python/muSpectre/tools.py +++ b/language_bindings/python/muSpectre/tools.py @@ -1,728 +1,718 @@ #!/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 from os import getcwd -import time def compute_real_space_lattice(resolutions, lengths, shift=np.array(None), d=np.array(None)): """ Set up the orthogonal real space lattice for a given resolution and cell length in each direction Parameters ---------- resolutions: np.array int Grid resolutions in the Cartesian directions. The returned grid have always the number of points given by resolution independent from 'd' For a three dimensional problem you need the following input form: [Nx, Ny, Nz] lengths: np.array float Physical size of the cell in the Cartesian directions. For a three dimensional problem you need the following input form: [Lx, Ly, Lz] shift: np.array float default: np.array(None), makes shift=np.zeros(ndim) All grid points are shifted by the array 'shift'. The shape has to be the spatial dimension (ndim) of the problem e.g 3D shift=np.array([1,2,3]). d: np.array float default np.arry(None), makes d=lengths/resolutions. Gives indirect the grid spacing by reducing the lengths by d, typically it is given by d=lengths/resolutions. For the default d, d is equal to the grid spacing. You should only change this option if it is really necessary! Return ------ X: np.array float Grid point position at each grid point """ res = np.array(resolutions) #grid resolution array lens = np.array(lengths) #box lengths array ndim = len(res) #dimension of the problem d = d #gridspacing in each direction #check defaults if shift.all() == None: shift = np.zeros(ndim) if d.all() == None: d = lens/res #set to normal gridspacing in each direction if ndim == 1: X = np.mgrid[0+shift[0]:lens[0]-d[0]+shift[0]:res[0]*1j].reshape(1,res[0]) elif ndim == 2: X = np.mgrid[0+shift[0]:lens[0]-d[0]+shift[0]:res[0]*1j, 0+shift[1]:lens[1]-d[1]+shift[1]:res[1]*1j ] elif ndim == 3: X = np.mgrid[0+shift[0]:lens[0]-d[0]+shift[0]:res[0]*1j, 0+shift[1]:lens[1]-d[1]+shift[1]:res[1]*1j, 0+shift[2]:lens[2]-d[2]+shift[2]:res[2]*1j ] return X def compute_reciprocal_space_lattice(resolutions, lengths, factor=2*np.pi): """ Set up the orthogonal reciprocal space lattice for a given resolution and cell length in each direction. Parameters ---------- resolutions: np.array int Grid resolutions in the Cartesian directions. For a three dimensional problem you need the following input form: [Nx, Ny, Nz] lengths: np.array float Physical size of the cell in the Cartesian directions. For a three dimensional problem you need the following input form: [Lx, Ly, Lz] factor: float default: 2*np.pi Factor with which the np.fft.fftfreq() frequencies are multiplied Return ------ q: np.array float Reciprocal vector at each grid point """ 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 q_1d = lambda i: factor*np.fft.fftfreq(res[i], d[i]) if ndim == 1: qx = q_1d(0) q = np.array(np.meshgrid(qx, indexing='ij')) elif ndim == 2: qx = q_1d(0) qy = q_1d(1) q = np.array(np.meshgrid(qx,qy, indexing='ij')) elif ndim == 3: qx = q_1d(0) qy = q_1d(1) qz = q_1d(2) q = np.array(np.meshgrid(qx, qy, qz, indexing='ij')) return q def extend_grid_by_pbc(grid_pbc, pbc_shift=np.array(None)): """ Extend a pbc grid in all dimensions by its periodic boundary (pb). The first entry in each direction is used to compute the new last entry. If the pb values differ by a shift from one side to the other, this shift can be given as pbc_shift. Parameters ---------- grid_pbc: np.ndarray float Grid, with periodic boundaries, which should be extended in each direction by the reconstructed periodic boundary. The first index 'i' is the value (value.shape=(dim)) of the grid at the grid coordinate 'x,y,z,...' with Nx, Ny, Nz,... grid points in each direction. Thus grid_pbc.shape = (dim,Nx,Ny,Nz,...). pbc_shift: np.array float default: np.array(None), makes shift=np.zeros(ndim) with ndim = grid.shape[0] The added periodic boundary values are shifted by the corresponding pbc_shift value. The shape of pbc_shift has to be the spatial dimension (ndim) of the problem e.g 3D pbc_shift=np.array([Lx,Ly,Lz]). Return ------ grid: np.array float Grid point position at each grid point """ ndim = grid_pbc.shape[0] #dimension of the problem #check default if pbc_shift.all() == None: pbc_shift = np.zeros(ndim) pad_w = ((0,0),)+((0,1),)*ndim grid = np.pad(grid_pbc, pad_width=pad_w, mode='wrap') #shift of the last value in each axis by pbc_shift in the corresponding #direction if ndim == 1: grid[0,-1] += pbc_shift[0] elif ndim == 2: grid[0,-1,:] += pbc_shift[0] grid[1,:,-1] += pbc_shift[1] elif ndim == 3: grid[0,-1,:,:] += pbc_shift[0] grid[1,:,-1,:] += pbc_shift[1] grid[2,:,:,-1] += pbc_shift[2] return grid def i_DFT(f_q, res, lens, x): """ Computes the inverse discrete Fourier transform (DFT) for given real space vectors x. + Internal the real space vector is separated into x = x_normal + delta_x. + There x_normal are the real space vectors one would expect in a normal ifft and + delta_x is the shift. Parameters ---------- f_q: np.ndarray, float or complex The Fourier coefficients (Fourier transformed of a real space function f). res: np.array, int Grid resolutions in the Cartesian directions. For a three dimensional problem you need the following input form: np.array([Nx, Ny, Nz]) lens: np.array, float Physical size of the cell in the Cartesian directions. For a three dimensional problem you need the following input form: np.array([Lx, Ly, Lz]) x: np.ndarray, float x vector used for the Fourier transformation with the shape (ndim, res). Return ------ f: np.ndarray, complex the inverse DFT of the signal f_q at the positions x. - - TODO - ---- - Seems to be very slow -> improve it... """ res = res lens = lens - ndim = len(res) + q = compute_reciprocal_space_lattice(res, lens) - f_q = f_q.reshape(f_q.shape + (1,)*ndim) + x_normal = compute_real_space_lattice(res, lens) + q_delta_x = (q * (x-x_normal)).sum(axis=0).reshape((1,)+tuple(res)) + f_q_shift = f_q * np.exp(1j * q_delta_x) - # qx (ndarray where each x vector is multiplied with the whole q on the grid) - shape = (1,) + 2*tuple(res) - qx = np.tensordot(q, x, axes=((0),(0))).reshape(shape) - ax = (0,) + tuple(np.arange(ndim+1,2*ndim+1)) + tuple(np.arange(1,ndim+1)) - f = 1/np.prod(res) * np.transpose(f_q * np.exp(1j*qx), axes=ax) \ - .sum( axis=tuple(-np.arange(1,ndim+1)) ) + f = np.fft.ifftn(f_q_shift, s=res) return f def compute_displacements(F, resolutions, lengths, order=0): """ Method to compute the real space displacement vectors 'displacements' 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] 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] 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 ------- initial_pos: np.ndarray, float initial positions X. displacements: np.ndarray, float displacement vectors pointing from the initial (undeformed) positions to the final (deformed) positions. Thus the final (deformed) structure is given by: initial_pos + displacements TODO ---- --> implement higher order corrections to an arbitrary order by the Lanczos sigma factor. """ 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 # Initialize the undeformed positions "X" and the q-vectors "q" X = compute_real_space_lattice(res, lens) q = compute_reciprocal_space_lattice(res, lens, factor=2*np.pi) ### 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 = F_homo.reshape((ndim,)*2+(1,)*ndim) F_nonhomo = F-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 # 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: # 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 #Keep in mind: This sets the integration constant to zero!!! u_q = - np.einsum("ij...,j...->i...", F_q, 1j*q)/q_norm_square elif order == 1: # first oder correction #reshape d as delta_x delta_x = d.reshape((ndim,)+(1,)*ndim) #compute the denominator and correct the zeros in it denom = (np.sin(q*delta_x)**2).sum(axis=0) denom.itemset(0,float("inf")) #correct for q=0 (Caution: This sets mean = 0!) #for even grids one also has to correct the sin for the gratest q wave length #because here you get sin(q_max)=0 ==> bad for the division! for i, n_grid_points in enumerate(res): if n_grid_points%2 == 0: #you have an even grid in direction i slice_pos = list((np.s_[:],)*ndim) slice_pos[i] = int(n_grid_points/2) slice_pos = tuple(slice_pos) denom[slice_pos] = float("inf") L_sigma = (1j*delta_x * np.sin(q*delta_x)) / denom u_q = - np.einsum("ij...,j...->i...", F_q, L_sigma) elif order == 2: # second oder correction #reshape d as delta_x delta_x = d.reshape((ndim,)+(1,)*ndim) sin_1 = 8*np.sin(q*delta_x) sin_1_square = (sin_1**2).sum(axis=0) sin_2 = np.sin(2*q*delta_x) sin_2_square = (sin_2**2).sum(axis=0) #compute the denominator and correct the zeros in it denom = (sin_1_square - sin_2_square) denom.itemset(0,float("inf")) #correct for q=0 (Caution: This sets mean = 0!) #for even grids one also has to correct the sin for the gratest q wave length #because here you get sin(q_max)=0 ==> bad for the division! for i, n_grid_points in enumerate(res): if n_grid_points%2 == 0: #you have an even grid in direction i slice_pos = list((np.s_[:],)*ndim) slice_pos[i] = int(n_grid_points/2) slice_pos = tuple(slice_pos) denom[slice_pos] = float("inf") L_sigma = (1j*6*delta_x*(sin_1+sin_2)) / denom u_q = - np.einsum("ij...,j...->i...", F_q, L_sigma) 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') 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) # only the real part is of interest # make the array C contiguous, without it is of mixed style nonaffine_displacements = np.ascontiguousarray( np.real(np.fft.ifftn(u_q, s=res)) ) initial_pos = X one = np.eye(ndim).reshape((ndim,)*2+(1,)*ndim) #'one tensor' affine_displacements = np.einsum("ij...,j...->i...", F_homo-one, X) displacements = affine_displacements + nonaffine_displacements return initial_pos, displacements def point_to_volume_grid(point_grid, resolutions, lengths, F, pbc = True, nonaffine_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 (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: np.ndarray float (with dimension of the problem) Grid around which the volume cells are constructed. The point_grid is assumed to be periodic, thus boundary points are occurring only once in the array. 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]. F: np.ndarray float The deformation gradient which lead to the deformed point_grid. pbc: bool default: True Decides if the returned volume_grid has periodic boundaries (surface points only on one side) (pbc=True) or non-periodic boundaries (pbc=False). nonaffine_displacements: np.ndarray float (with dimension of the problem) default: None (The nonaffine displacements are reconstructed from the deformation gradient F and the point_grid.) 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: np.ndarray float (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 ---- --> make an option to decide whether the deformation is handled periodically (a fixed box and points which flow out coming in again on the periodic boundary) or the deformation is shown as in real space (deformation can go out of the box). --> maybe introduce an default for 'F' e.g. np.eye() ? (If one does not know the deformation gradient or the deformation was zero...) """ #initialization res = np.array(resolutions) lens = np.array(lengths) d = lens/res ndim = len(res) #homogeneous deformation gradient, F_homo 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) #non-affine deformation, nonaff_d nonaff_d = nonaffine_displacements #check default if nonaff_d.all() == None: #construct nonaff_d from: nonaff_d = point_grid - F_homo*X_undeformed X_und = compute_real_space_lattice(res, lens) #undeformed (und) grid nonaff_d = point_grid - np.einsum('ij...,j...->i...', F_homo, X_und) nonaff_d_sum = np.zeros(point_grid.shape) weight = 0 ### weighted sum of nonaff_d; # The position of an edge is given as the homogeneous deformation (X_vol= # F_homo*X_shift) of the undeformed edge grid (X_shift), plus the average # nonaffine deformation of the edge point surrounding grid points. for dim in range(ndim+1): ## iterate over the dimensions of the problem for ax in itertools.combinations(range(1,ndim+1),dim): ## iterate over all neigbors of an volume edge point ## 8 in 3D, 4 in 2D, 2 in 1D nonaff_d_sum += np.roll(nonaff_d, 1, axis=ax) weight += 1 #average nonaffine deformation of an edge point nonaff_d_av_pbc = 1/weight * nonaff_d_sum if pbc: #pbc=True, return a periodic boundary volume grid #(surface points only on one side) #edge positions in the undeformed grid X_shift_pbc = compute_real_space_lattice(res, lens, shift=-d/2) #X_vol, edge positions in the homogeneous deformed grid X_vol_pbc = np.einsum('ij...,j...->i...',F_homo, X_shift_pbc) volume_grid = X_vol_pbc + nonaff_d_av_pbc - #return volume_grid_pbc - else: #pbc=False, return a grid with full boundaries on all sides #edge positions in the undeformed grid X_shift = compute_real_space_lattice(res+1, lens+d, shift=-d/2, d=d) #X_vol, edge positions in the homogeneous deformed grid X_vol = np.einsum('ij...,j...->i...',F_homo, X_shift) nonaff_d_av = extend_grid_by_pbc(nonaff_d_av_pbc) volume_grid = X_vol + nonaff_d_av return volume_grid def point_to_volume_grid_fourier(point_grid, resolutions, lengths, F, pbc = True): """ Function which converts the point grid (grid on which computations are done) into a volume grid (for visualization). The interpolation is done by Fourier interpolate the grid onto the volume edge points. The grid is assumed to be periodic. Parameters ---------- point_grid: np.ndarray float (with dimension of the problem) Grid around which the volume cells are constructed. The point_grid is assumed to be periodic, thus boundary points are occurring only once in the array. 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]. F: np.ndarray float The deformation gradient which lead to the deformed point_grid. pbc: bool default: True Decides if the returned volume_grid has periodic boundaries (surface points only on one side) (pbc=True) or non-periodic boundaries (pbc=False). Returns ------- volume_grid: np.ndarray float (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 ---- --> make an option to decide whether the deformation is handled periodically (a fixed box and points which flow out coming in again on the periodic boundary) or the deformation is shown as in real space (deformation can go out of the box). """ #initialization res = np.array(resolutions) lens = np.array(lengths) d = lens/res ndim = len(res) ### Fourier interpolate the non-affine deformation: #Interpolation grid x_interpol = compute_real_space_lattice(resolutions = res, lengths = lens, shift = -d/2, d = d ) x_initial = compute_real_space_lattice(resolutions = res, lengths = lens) #Compute affine and non-affine deformation dim_axis = tuple(-np.arange(1,ndim+1)) #the last ndim entries F_homo = 1./(np.prod(res)) * F.sum(axis=dim_axis) affine_d_initial = np.einsum("ij...,j...->i...", F_homo, x_initial) affine_d_interpol = np.einsum("ij...,j...->i...", F_homo, x_interpol) non_affine_d = point_grid - affine_d_initial #Fourier coefficients of non-affine deformation f_q = np.fft.fftn(non_affine_d, s=res) #Fourier interpolation - st = time.time() volume_grid = np.real(i_DFT(f_q, res, lens, x_interpol)) + affine_d_interpol - et = time.time() - print('time for Fourier trafo: ', et-st) if not pbc: #Non periodic boundaries, map the boundary points on both sides of the array #affine_d = extend_grid_by_pbc(affine_d, pbc_shift=lens) volume_grid = extend_grid_by_pbc(volume_grid, pbc_shift=lens) return volume_grid def write_structure_as_vtk(file_name, positions, cellData=None, pointData=None): """ Function to save a structure as vts file by using pyevtk.hl. The file: file_name.vts is written. You can open it by e.g. the software Paraview. Parameters ---------- file_name: str name of the 'vts' file which will be written. positions: np.ndarray, float grid point positions with the shape (dim, Nx, Ny, Nz). Where 'dim' is the dimension (1D, 2D or 3D) and Nx,Ny,Nz are the number of grid points in each direction. cellData: dictionary default: None, np.zeros of the right shape will be constructed. data associated with voxel properties. You can give a dictionary of arbitrary length. The grid shape of the entries should correspond to (Nx-1, Ny-1, Nz-1). pointData: dictionary default: None, np.zeros of the right shape will be constructed. data associated to each point / point properties. You can give a dictionary of arbitrary length. The grid shape of the entries should correspond to (Nx, Ny, Nz). Returns ------- Writes a '.vts' file with the file name, file_name.vts. If it ends successfully it returns the stirng 'The nD grid: file_name.vts was written.' """ from pyevtk.hl import gridToVTK #initialize parameters dim = len(positions.shape)-1 cdata = cellData pdata = pointData grid = np.array(positions.shape[1:]) #number of grid points in each direction #set positions and reshape pdata and cdata if necessary if dim == 1: x = np.ascontiguousarray(positions[0].reshape((len(positions[0]),1,1))) y = np.zeros(x.shape) z = np.zeros(x.shape) init_d = 1 #initial dimension of the data vol_shape = (grid[0]-1, 1, 1) #expected shape for cell data elif dim == 2: x = np.expand_dims(positions[0],axis=2) y = np.expand_dims(positions[1],axis=2) z = np.zeros(x.shape) init_d = 2 vol_shape = (grid[0]-1, grid[1]-1, 1) elif dim == 3: x = positions[0] y = positions[1] z = positions[2] init_d = 3 vol_shape = tuple(grid - 1) #check pdata and cdata if arrays are 3D, because gridToVTK expects 3D data if (cdata != None): for k in cdata.keys(): if len(cdata[k].shape) >= 3: if cdata[k].shape[-3:] == x.shape: print('\n\nYour cell data has the shape of your positions {}.\n' 'This means you loose the last entry in each direction.\n' 'You should better have the shape {}.' .format(x.shape,vol_shape)) elif cdata[k].shape[-3:] != vol_shape: print('\n\nThere is probably a problem:\n'\ 'cell data "{}" shape is {} and your grid is {}.' .format(k, cdata[k].shape, x.shape)) new_shape = cdata[k].shape + (1,)*(3-init_d) cdata[k] = cdata[k].reshape(new_shape) print('I reshaped the data to {}, to write the vts file.' .format(new_shape)) else: print('\n\nI try to fix the entry "{}" with shape {} in your cell data.' .format(k, cdata[k].shape)) new_shape = cdata[k].shape + (1,)*(3-len(cdata[k].shape)) cdata[k] = cdata[k].reshape(new_shape) print('The data was reshaped to {}'.format(new_shape)) #check if cell data has now the right shape for i, n in enumerate(cdata[k].shape[-3:]): if n != vol_shape[i]: print('\n\nWARNING:\nYour cell data "{}" with shape {} seem to be ' 'wrong.\nI would have expected a grid (last entries in cell ' 'data shape) with shape {},\ncorresponding to your input grid ' 'with shape {}.\nPlease have a look!' .format(k, cdata[k].shape, vol_shape, x.shape)) break if (pdata != None): for k in pdata.keys(): if len(pdata[k].shape) >= 3: if pdata[k].shape[-3:] != x.shape: print('\n\nThere is probably a problem:\n'\ 'point data "{}" shape is {} and your grid is {}.' .format(k, pdata[k].shape, x.shape)) new_shape = pdata[k].shape + (1,)*(3-init_d) pdata[k] = pdata[k].reshape(new_shape) print('I reshaped the data to {}, to write the vts file.' .format(new_shape)) else: print('\n\nI try to fix the entry "{}" with shape {} in your point data.' .format(k, pdata[k].shape)) new_shape = pdata[k].shape + (1,)*(3-len(pdata[k].shape)) pdata[k] = pdata[k].reshape(new_shape) print('The data was reshaped to {}'.format(new_shape)) #check if point data has now the right shape for i, n in enumerate(pdata[k].shape[-3:]): if n != x.shape[i]: print('\n\nWARNING:\nYour point data "{}" with shape {} seem to be ' 'wrong.\nI would have expected a grid (last entries in point ' 'data shape) with shape {},\ncorresponding to your input grid ' 'with shape {}.\nPlease have a look!' .format(k, cdata[k].shape, x.shape, x.shape)) break #default values if cdata == None: cdata = {'no data available' : np.zeros(x.shape)} if pdata == None: pdata = {'no data available' : np.zeros(x.shape)} gridToVTK(file_name, x, y, z, cellData=cdata, pointData=pdata) cwd = getcwd()+'/' message = '\nThe {}D grid:\n{}_grid.vts\nwas written.\n'.format(dim, cwd+file_name) return message diff --git a/tests/python_muSpectre_tools_test.py b/tests/python_muSpectre_tools_test.py index cff8712..44d83e6 100644 --- a/tests/python_muSpectre_tools_test.py +++ b/tests/python_muSpectre_tools_test.py @@ -1,750 +1,753 @@ #!/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 scipy.misc as sm -import itertools import tempfile import os -import time 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)! 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 and 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 ---- * tests for: -> compute_real_space_lattice() -> compute_reciprocal_space_lattice() -> extend_grid_by_pbc() * better tests for: -> test_point_to_volume_grid() (off diagonal F and nonaff displacements) """ def setUp(self): self.lengths = np.array([2.4, 3.7, 4.1]) - self.resolutions = np.array([4, 19, 13]) + self.resolutions = np.array([32, 19, 51]) #tolerance for the norm(Chi_X - integral(F)) self.norm_tol = 1e-8 #tollerance for the central difference approximation test self.norm_tol_cent_diff = 1e-8 #tollerance for i_DFT test self.norm_tol_i_DFT = 1e-8 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)=sin(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_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) + X = µ.tools.compute_real_space_lattice(res, lens) - #real signal + ### 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 = np.cos(2*np.pi/lens.reshape((ndim,)+(1,)*ndim) * X) + ### 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_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 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. """ ###### test for the deformation gradient integration scheme ###### ### F=1I, 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] initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### cosine diagonal deformation gradient # The resolution 30 and legth 2.2 in x-direction test a special case of # grid, here X = np.mgrid[0:lens[0]:d[0], 0:lens[1]:d[1], 0:lens[2]:d[2]] # fails because it makes 31 grid points in x direction. res = [4, 30, 2] lens = [4, 2.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]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (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 = [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 initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), 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 # #initial_pos, displacements = µ.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) F[0,0,:] = 5 + 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*X[0]) Chi_X[0,:] = 5*X[0] + 2*np.sin(2*np.pi*X[0]/lens[0]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### 2D ### res = [5, 5] lens = [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[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]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### 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]) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=0) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ###### tests for higher order corrections ###### order_o = [1, 2] ### cosine diagonal deformation gradient - res = [6, 10] + res = [30, 30] 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,:,:] = 0.2*2*np.pi/lens[i] * np.cos(2*np.pi*X[i]/lens[i]) - Chi_X[i,:,:] = 0.2 * np.sin(2*np.pi*X[i]/lens[i]) + 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]) # 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) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (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 = [2000, 1, 1] 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,:,:,:] = 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 # 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) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### 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,:,:] = 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]) # 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) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (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,:,:,:] = 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]) # 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) initial_pos, displacements = µ.tools.compute_displacements(F, res, lens, order=n) self.assertLess(np.linalg.norm(Chi_X - (initial_pos + displacements)), self.norm_tol) ### real test, shear of material with two different hardnesses/ Young moduli. ### Can I see Gibbs ringing? #initialize res = [101, 101, 1] lens = [10, 10, 1] ndim = len(res) #dimension of the problem Nx, Ny, Nz = res Lx, Ly, Lz = 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 #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,Nz), dtype=int) #0 for surrounding matrix boundary = res[1]//2 phase[:, boundary:, :] = 1 cell = µ.Cell(res, lens, formulation) mat = µ.material.MaterialLinearElastic4_3d.make(cell, "material") for i, pixel in enumerate(cell): #add Young and Poisson depending on the material index m_i = phase.flatten()[i] #m_i = material index (phase index) mat.add_pixel(pixel, Young[m_i], Poisson[m_i]) cell.initialise() DelF = np.array([[0 , 0.1, 0], [0.1, 0 , 0], [0 , 0 , 0]]) solver_newton = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose) r_newton = µ.solvers.newton_cg(cell, DelF, solver_newton, newton_tol, equil_tol, verbose) ### check strains strain = r_newton.grad.T.reshape(*res, 3, 3).transpose((4,3,0,1,2)) one = np.eye(ndim).reshape((ndim,)*2+(1,)*ndim) F = strain + one initial_pos, displacements = \ µ.tools.compute_displacements(F, res, lens, order) fin_pos = initial_pos + displacements ### analytic solution # Why L - gs / 2 ?! # The analytic signal doesn't know about the periodic boundary conditions, # 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) #(Ly - gs[1]) / 2 , analytic height l_hard_ana = gs[1]*(Ny - boundary - 1/2) #(Ly - gs[1]) / 2 , analytic height l_soft = gs[1] * boundary #height material soft l_hard = gs[1] * (Ny-boundary) #height material hard mu = np.array(Young) / (2 * (1+np.array(Poisson))) #shear modul (mu_soft, mu_hard) 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)#X.copy() #initialize new positions #soft d[0, :, :boundary, :] += g_soft/2 * X[1, :, :boundary, :] #d_soft x coordinate d[1, :, :boundary, :] += gamma/2 * X[0, :, :boundary, :] #d_soft y coordinate #hard d[0, :, boundary:, :] += g_hard/2 * (X[1, :, boundary:, :]-l_soft_ana) \ + g_soft/2 * l_soft_ana d[1, :, boundary:, :] += gamma/2 * (X[0, :, boundary:, :]) #correct by integration constant # shift the analytic solution such that: # -> the average nonaffine deformation is zero (integral of the nonaffine # deformation gradient + integration const). # ==> choose the right 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,0] / Ny #final positions = displacements + old positions + integration const. d += X + np.array([0, int_const, 0]).reshape((ndim,)+(1,)*ndim) #self.assertLess(np.linalg.norm(fin_pos - d), self.norm_tol) ### real test, integrate a delta function deformation gradient to a step function. ### Can I see Gibbs ringing? # does this test make sense? -> discuss with Till! 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???!!! """ - st = time.time() - 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 ) - et = time.time() - print('By hand interpol {}: {} s'.format(res[:dim], et-st)) - 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(). """ - st = time.time() - 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])) #print('dimension: ', dim) #print('res: ', 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 ) - et = time.time() - print('Fourier interpol {}: {} s'.format(res[:dim], et-st)) - 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()