diff --git a/bin/visualization_example.py b/bin/visualization_example.py new file mode 100644 index 0000000..ebef427 --- /dev/null +++ b/bin/visualization_example.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +""" +@file visualization_example.py + +@author Richard Leute + +@date 14 Jan 2019 + +@brief small example of how one can use the visualization tools: + gradient_integration() and vtk_export() + +@section LICENSE + +Copyright © 2019 Till Junge, Richard Leute + +µ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 µSpectre; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 59 Temple Place - Suite 330, +Boston, MA 02111-1307, USA. + +Additional permission under GNU GPL version 3 section 7 + +If you modify this Program, or any covered work, by linking or combining it +with proprietary FFT implementations or numerical libraries, containing parts +covered by the terms of those libraries' licenses, the licensors of this +Program grant you additional permission to convey the resulting work. +""" + + +import numpy as np +import sys +import os +sys.path.append(os.path.join(os.getcwd(), "../build/language_bindings/python/")) +import muSpectre as µ +import muSpectre.gradient_integration as gi +import muSpectre.vtk_export as vt_ex + + +### Input parameters ### +#----------------------# +#general +resolutions = [71, 71, 3] +lengths = [1.0, 1.0, 0.2] +Nx, Ny, Nz = resolutions +formulation = µ.Formulation.finite_strain +Young = [10, 20] #Youngs modulus for each phase +Poisson = [0.3, 0.4] #Poissons ratio for each phase + +#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 + +# sinusoidal bump +d = 10 #thickness of the phase with higher viscosity in pixles +l = Nx//3 #length of bump +h = Ny//4 #height of bump + +low_y = (Ny-d)//2 #lower y-boundary of phase +high_y = low_y+d #upper y-boundary of phase + +left_x = (Nx-l)//2 #boundaries in x direction left +right_x = left_x + l #boundaries in x direction right +x = np.arange(l) +p_y = h*np.sin(np.pi/(l-1)*x) #bump function + +xy_bump = np.ones((l,h,Nz)) #grid representing bumpx +for i,treshold in enumerate(np.round(p_y)): + xy_bump[i,int(treshold):,:] = 0 + +phase = np.zeros(resolutions, dtype=int) #0 for surrounding matrix +phase[:, low_y:high_y, :] = 1 +phase[left_x:right_x, high_y:high_y+h, :] = xy_bump + + +### Run muSpectre ### +#-------------------# +cell = µ.Cell(resolutions, + lengths, + 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() #initialization of fft to make faster fft +DelF = np.array([[0 , 0.7, 0], + [0 , 0 , 0], + [0 , 0 , 0]]) + +solver_newton = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose) +result = µ.solvers.newton_cg(cell, DelF, solver_newton, + newton_tol, equil_tol, verbose) + +# print solver results +print('\n\nstatus messages of OptimizeResults:') +print('-----------------------------------') +print('success: ', result.success) +print('status: ', result.status) +print('message: ', result.message) +print('# iterations: ', result.nb_it) +print('# cell evaluations: ', result.nb_fev) +print('formulation: ', result.formulation) + + +### Visualization ### +#-------------------# +#integration of the deformation gradient field +placement_n, x = gi.compute_placement(result, lengths, resolutions, order = 0) + +### some fields which can be added to the visualization +#2-tensor field containing the first Piola Kirchhoff stress +PK1 = gi.reshape_gradient(result.stress, resolutions) +#scalar field containing the distance to the origin O +distance_O = np.linalg.norm(placement_n, axis=-1) + +#random fields +center_shape = tuple(np.array(x.shape[:-1])-1)#shape of the center point grid +dim = len(resolutions) +scalar_c = np.random.random(center_shape) +vector_c = np.random.random(center_shape + (dim,)) +tensor2_c = np.random.random(center_shape + (dim,dim)) +scalar_n = np.random.random(x.shape[:-1]) +vector_n = np.random.random(x.shape[:-1] + (dim,)) +tensor2_n = np.random.random(x.shape[:-1] + (dim,dim)) + +#write dictionaries with cell data and point data +c_data = {"scalar field" : scalar_c, + "vector field" : vector_c, + "2-tensor field" : tensor2_c, + "PK1 stress" : PK1, + "phase" : phase} +p_data = {"scalar field" : scalar_n, + "vector field" : vector_n, + "2-tensor field" : tensor2_n, + "distance_O" : distance_O} + +vt_ex.vtk_export(fpath = "visualization_example", + x_n = x, + placement = placement_n, + point_data = p_data, + cell_data = c_data) diff --git a/language_bindings/python/muSpectre/gradient_integration.py b/language_bindings/python/muSpectre/gradient_integration.py index bca673d..528a5d6 100644 --- a/language_bindings/python/muSpectre/gradient_integration.py +++ b/language_bindings/python/muSpectre/gradient_integration.py @@ -1,366 +1,365 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file gradient_integration.py @author Till Junge @date 22 Nov 2018 @brief Functions for the integration of periodic first- and second-rank tensor fields on an n-dimensional rectangular grid Copyright © 2018 Till Junge µSpectre is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser 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 Lesser General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ import numpy as np import sys def compute_wave_vectors(lengths, resolutions): """Computes the wave vectors for a dim-dimensional rectangular or cubic (or hypercubic) grid as a function of the edge lengths and resolutions. Note: the norm of the wave vectors corresponds to the angular velocity, not the frequency. Keyword Arguments: lengths -- np.ndarray of length dim with the edge lengths in each spatial dimension (dtype = float) resolutions -- np.ndarray of length dim with the resolutions in each spatial dimension (dtype = int) Returns: np.ndarary of shape resolutions + [dim]. The wave vector for a given pixel/voxel is given in the last dimension """ return np.moveaxis(np.meshgrid( *[2*np.pi*np.fft.fftfreq(r, l/r) for l,r in zip(lengths, resolutions)], indexing="ij"), 0, -1) def compute_grid(lengths, resolutions): """For a dim-dimensional pixel/voxel grid, computes the pixel/voxel centre and corner positions as a function of the grid's edge lengths and resolutions Keyword Arguments: lengths -- np.ndarray of length dim with the edge lengths in each spatial dimension (dtype = float) resolutions -- np.ndarray of length dim with the resolutions in each spatial dimension (dtype = int) Returns: tuple((x_n, x_c)) two ndarrays with nodal/corner positions and centre positions respectively. x_n has one more entry in every direction than the resolution of the grid (added points correspond to the periodic repetitions) """ x_n = np.moveaxis(np.meshgrid( *[np.linspace(0, l, r+1) for l,r in zip(lengths, resolutions)], indexing="ij"), 0 ,-1) dx = lengths/resolutions x_c = np.moveaxis(np.meshgrid( *[np.linspace(0, l, r, endpoint=False) for l,r in zip(lengths, resolutions)], indexing="ij"), 0, -1) + .5*dx return x_n, x_c def reshape_gradient(F, resolutions): """reshapes a flattened second-rank tensor into a multidimensional array of shape resolutions + [dim, dim]. Note: this reshape entails a copy, because of the column-major to row-major transposition between Eigen and numpy Keyword Arguments: F -- flattenen array of gradients as in OptimizeResult resolutions -- np.ndarray of length dim with the resolutions in each spatial dimension (dtype = int) Returns: np.ndarray """ dim = len(resolutions) if not isinstance(resolutions, list): raise Exception("resolutions needs to be in list form, "+ "for concatenation") expected_input_shape = [np.prod(resolutions) * dim**2] output_shape = resolutions + [dim, dim] if not ((F.shape[0] == expected_input_shape[0]) and (F.size == expected_input_shape[0])): raise Exception("expected gradient of shape {}, got {}".format( expected_input_shape, F.shape)) order = list(range(dim+2)) order[-2:] = reversed(order[-2:]) return F.reshape(output_shape).transpose(*order) def complement_periodically(array, dim): """Takes an arbitrary multidimensional array of at least dimension dim and returns an augmented copy with periodic copies of the left/lower entries in the added right/upper boundaries Keyword Arguments: array -- arbitrary np.ndarray of at least dim dimensions dim -- nb of dimension to complement periodically """ shape = list(array.shape) shape[:dim] = [d+1 for d in shape[:dim]] out_arr = np.empty(shape, dtype = array.dtype) sl = tuple([slice(0, s) for s in array.shape]) out_arr[sl] = array for i in range(dim): lower_slice = tuple([slice(0,s) if (d != i) else 0 for (d,s) in enumerate(shape)]) upper_slice = tuple([slice(0,s) if (d != i) else -1 for (d,s) in enumerate(shape)]) out_arr[upper_slice] = out_arr[lower_slice] return out_arr def get_integrator(x, freqs, order=0): """returns the discrete Fourier-space integration operator as a function of the position grid (used to determine the spatial dimension and resolution), the wave vectors, and the integration order Keyword Arguments: x -- np.ndarray of pixel/voxel centre positons in shape resolution + [dim] freqs -- wave vectors as computed by compute_wave_vectors order -- (default 0) integration order. 0 stands for exact integration Returns: (dim, shape, integrator) """ dim = x.shape[-1] shape = x.shape[:-1] delta_x = x[tuple([1]*dim)] - x[tuple([0]*dim)] def correct_denom(denom): """Corrects the denominator to avoid division by zero for freqs = 0 and on even grids for freqs*delta_x=pi.""" denom[tuple([0]*dim)] = 1. for i, n in enumerate(shape): if n%2 == 0: denom[tuple([np.s_[:]]*i + [n//2] + [np.s_[:]]*(dim-1-i))] = 1. return denom[..., np.newaxis] if order == 0: freqs_norm_square = np.einsum("...i, ...i -> ...", freqs, freqs) freqs_norm_square.reshape(-1)[0] = 1. integrator = 1j * freqs/freqs_norm_square[...,np.newaxis] # Higher order corrections after: # A. Vidyasagar et al., Computer Methods in Applied Mechanics and # Engineering 106 (2017) 133-151, sec. 3.4 and Appendix C. # and # P. Eisenlohr et al., International Journal of Plasticity 46 (2013) # 37-53, Appendix B, eq. 23. elif order == 1: sin_1 = np.sin(freqs*delta_x) denom = correct_denom((sin_1**2).sum(axis=-1)) integrator = 1j*delta_x*sin_1 / denom elif order == 2: sin_1, sin_2 = 8*np.sin(freqs*delta_x), np.sin(2*freqs*delta_x) denom = correct_denom((sin_1**2).sum(axis=-1) - (sin_2**2).sum(axis=-1)) integrator = 1j*6*delta_x*(sin_1+sin_2) / denom else: raise Exception("order '{}' is not implemented".format(order)) return dim, shape, integrator def integrate_tensor_2(grad, x, freqs, staggered_grid=False, order=0): """Integrates a second-rank tensor gradient field to a chosen order as a function of the given field, the grid positions, and wave vectors. Optionally, the integration can be performed on the pixel/voxel corners (staggered grid). Keyword Arguments: grad -- np.ndarray of shape resolution + [dim, dim] containing the second-rank gradient to be integrated x -- np.ndarray of shape resolution + [dim] (or augmented resolution + [dim]) containing the pixel/voxel centre positions (for un-staggered grid integration) or the pixel /voxel corner positions (for staggered grid integration) freqs -- wave vectors as computed by compute_wave_vectors staggered_grid -- (default False) if set to True, the integration is performed on the pixel/voxel corners, rather than the centres. This leads to a different integration scheme order -- (default 0) integration order. 0 stands for exact integration Returns: np.ndarray contaning the integrated field """ dim, shape, integrator = get_integrator(x, freqs, order) axes = range(dim) grad_k = np.fft.fftn(grad, axes=axes) f_k = np.einsum("...j, ...ij -> ...i", integrator, grad_k) normalisation = np.prod(grad.shape[:dim]) grad_k_0 = grad_k[tuple((0 for _ in range(dim)))].real/normalisation homogeneous = np.einsum("ij, ...j -> ...i", grad_k_0, x) if not staggered_grid: fluctuation = -np.fft.ifftn(f_k, axes=axes).real else: del_x = (x[tuple((1 for _ in range(dim)))] - x[tuple((0 for _ in range(dim)))]) k_del_x = np.einsum("...i, i ->...", freqs, del_x)[...,np.newaxis] if order == 0: shift = np.exp(-1j * k_del_x/2) elif order == 1: shift = (np.exp(-1j * k_del_x) + 1) / 2 elif order == 2: shift = np.exp(-1j*k_del_x/2) * np.cos(k_del_x/2) *\ (np.cos(k_del_x) - 4) / (np.cos(k_del_x/2) - 4) fluctuation = complement_periodically( -np.fft.ifftn(shift*f_k, axes=axes).real, dim) return fluctuation + homogeneous def integrate_vector(df, x, freqs, staggered_grid=False, order=0): """Integrates a first-rank tensor gradient field to a chosen order as a function of the given field, the grid positions, and wave vectors. Optionally, the integration can be performed on the pixel/voxel corners (staggered_grid) Keyword Arguments: df -- np.ndarray of shape resolution + [dim] containing the first-rank tensor gradient to be integrated x -- np.ndarray of shape resolution + [dim] (or augmented resolution + [dim]) containing the pixel/voxel centre positions (for un-staggered grid integration) or the pixel /voxel corner positions (for staggered grid integration) freqs -- wave vectors as computed by compute_wave_vectors staggered_grid -- (default False) if set to True, the integration is performed on the pixel/voxel corners, rather than the centres. This leads to a different integration scheme... order -- (default 0) integration order. 0 stands for exact integration Returns: np.ndarray contaning the integrated field """ dim, shape, integrator = get_integrator(x, freqs, order) axes = range(dim) df_k = np.fft.fftn(df, axes=axes) f_k = np.einsum("...i, ...i -> ...", df_k, integrator) df_k_0 = df_k[tuple((0 for _ in range(dim)))].real homogeneous = np.einsum("i, ...i -> ...", df_k_0, x/np.prod(shape)) if not staggered_grid: fluctuation = -np.fft.ifftn(f_k, axes=axes).real else: del_x = x[tuple((1 for _ in range(dim)))] \ - x[tuple((0 for _ in range(dim)))] k_del_x = np.einsum("...i, i ->...", freqs, del_x) if order == 0: shift = np.exp(-1j * k_del_x/2) elif order == 1: shift = (np.exp(-1j * k_del_x) + 1) / 2 elif order == 2: shift = np.exp(-1j*k_del_x/2) * np.cos(k_del_x/2) *\ (np.cos(k_del_x) - 4) / (np.cos(k_del_x/2) - 4) fluctuation = complement_periodically( -np.fft.ifftn(shift*f_k, axes=axes).real, dim) return fluctuation + homogeneous def compute_placement(result, lengths, resolutions, order=0, formulation=None): """computes the placement (the sum of original position and displacement) as a function of a OptimizeResult, domain edge lengths, domain discretisation resolutions, the chosen integration order and the continuum mechanics description(small or finite strain description) Keyword Arguments: result -- OptimiseResult, or just the grad field of an OptimizeResult lengths -- np.ndarray of length dim with the edge lengths in each spatial dimension (dtype = float) resolutions -- np.ndarray of length dim with the resolutions in each spatial dimension (dtype = int) order -- (default 0) integration order. 0 stands for exact integration formulation -- (default None) the formulation is derived from the OptimiseResult argument. If this is not possible you have to fix the formulation to either "small_strain" or "finite_strain". (dtype = str) Returns: (placement, x_n) tuple of ndarrays containing the placement and the corresponding original positions """ lengths = np.array(lengths) resolutions = np.array(resolutions) - dim = resolutions.shape[0] #Check whether result is a np.array or an OptimiseResult object if isinstance(result, np.ndarray): if formulation == None: #exit the program, if the formulation is unknown! raise ValueError('\n' 'You have to specify your continuum mechanics description.\n' 'Either you use a formulation="small_strain" or ' '"finite_strain" description.\n' 'Otherwise you can give a result=OptimiseResult object, which ' 'tells me the formulation.') form = formulation grad = reshape_gradient(result, resolutions.tolist()) else: form = str(result.formulation)[12:] if form != formulation and formulation != None: #exit the program, if the formulation is ambiguous! raise ValueError('\nThe given formulation "{}" differs from the ' 'one saved in your result "{}"!' .format(formulation, form)) grad = reshape_gradient(result.grad, resolutions.tolist()) #reshape the gradient depending on the formulation if form == "small_strain" or form == "small_strain_sym": raise NotImplementedError('\nIntegration of small strains' 'is not implemented yet!') elif form == "finite_strain": - grad = grad - np.eye(dim).reshape((1,)*dim + (dim,)*2) + grad = grad else: raise ValueError('\nThe formulation: "{}" is unknown!' .format(formulation)) #compute the placement by integrating x_n, x_c = compute_grid(lengths, resolutions) freqs = compute_wave_vectors(lengths, resolutions) placement = integrate_tensor_2(grad, x_n, freqs, staggered_grid=True, order=order) return placement, x_n diff --git a/tests/python_muSpectre_gradient_integration_test.py b/tests/python_muSpectre_gradient_integration_test.py index d74d4ca..7fcc399 100644 --- a/tests/python_muSpectre_gradient_integration_test.py +++ b/tests/python_muSpectre_gradient_integration_test.py @@ -1,661 +1,661 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file python_muSpectre_gradient_integration_test.py @author Richard Leute @date 23 Nov 2018 @brief test the functionality of gradient_integration.py @section LICENSE Copyright © 2018 Till Junge, Richard Leute µ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 µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ import unittest import numpy as np import scipy.misc as sm import itertools from python_test_imports import µ ### Helper functions def init_X_F_Chi(lens, res, rank=2): """ Setup all the needed parameters for initialization of the deformation gradient F and the corresponding deformation map/field Chi_X. Keyword Arguments: lens -- list [Lx, Ly, ...] of box lengths in each direction (dtype=float) res -- list [Nx, Ny, ...] of grid resoultions (dtype = int) rank -- int (default=2), rank of the deformation gradient tensor F. (dtype = int) Returns: d : np.array of grid spacing for each spatial direction (dtype = float) dim : int dimension of the structure, derived from len(res). x_n : np.ndarray shape=(res.shape+1, dim) initial nodal/corner positions as created by gradient_integration.compute_grid (dtype = float) x_c : np.ndarray shape=(res.shape+1, dim) initial cell center positions as created by gradient_integration.compute_grid (dtype = float) F : np.zeros shape=(res.shape, dim*rank) initialise deformation gradient (dtype = float) Chi_n: np.zeros shape=((res+1).shape, dim) initialise deformation field (dtype = float) Chi_c: np.zeros shape=(res.shape, dim) initialise deformation field (dtype = float) freqs: np.ndarray as returned by compute_wave_vectors(). (dtype = float) """ lens = np.array(lens) res = np.array(res) d = lens / res dim = len(res) x_n, x_c = µ.gradient_integration.compute_grid(lens, res) F = np.zeros(x_c.shape + (dim,)*(rank-1)) Chi_n = np.zeros(x_n.shape) Chi_c = np.zeros(x_c.shape) freqs = µ.gradient_integration.compute_wave_vectors(lens, res) return d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs def correct_to_zero_mean(Chi, d, nodal=False): """ corrects the displacement field such that it's integral is zero. By this one can get rid of a constant factor in the deformation gradient. This function is specialized for this file and should not be used somewhere else. Keywords: Chi : np.ndarray of the uncorrected analytic placements (dtype = float) d : np.array of the gridspacing in each spatial direction (dtype = float) nodal : bool (default False) specifies if the input are nodal or cell/center values. Default interpretation are cell/center values. Returns: Chi : np.ndarray of zero mean corrected analytic placements (dtype = float) """ Chi_zm = np.copy(Chi) res = np.array(Chi_zm.shape[:-1]) dim = res.size Chi_zm -= (Chi_zm.sum(axis=tuple(range(dim)))/np.prod(res))\ .reshape((1,)*dim + (dim,)) return Chi_zm def test_integrate(order, F, Chi_n, Chi_c, tol): """ make the convergence tests for the integration Keywords: order : list of integration orders which are tested (dtype = int) F : np.ndarray applied deformation gradient (dtype = float) Chi_n : np.ndarray expected node positions (dtype = float) Chi_c : np.ndarray expected cell positions (dtype = float) tol : list of tolerances for each order. If it is a single value the same tolerance is used for each order. (dtype = float) """ print('Maybe implement a function like this...') def central_diff_derivative(data, d, order, rank=1): """ Compute the first derivative of a function with values 'data' with the central difference approximation to the order 'order'. The function values are on a rectangualar grid with constant grid spacing 'd' in each direction. CAUTION: The function is assumed to be periodic (pbc)! Thus, if there is a discontinuity at the boundaries you have to expect errors in the derivative at the vicinity close to the discontinuity. Keyword Arguments: data -- np.ndarray of shape=(resolution, dim*rank) function values on an equally spaced grid, with grid spacing 'd' (dtype = float) d -- scalar or np.array of grid spacing in each direction. Scalar is interpreted as equal spacing in each direction (dtype = float) order -- int >= 1, gives the accuracy order of the central difference approximation (dtype = int) rank -- int, rank of the data tensor Returns: deriv: np.ndarray of shape=(resolution, dim, dim) central difference derivative of given order (dtype = float) """ dim = len(data.shape)-rank weights = sm.central_diff_weights(2*order + 1) deriv = np.zeros(data.shape + (dim,)) for ax in range(dim): for i in range(2*order + 1): deriv[...,ax] += weights[i]*np.roll(data, order-i, axis=ax) / d[ax] return deriv class MuSpectre_gradient_integration_Check(unittest.TestCase): """ Check the implementation of all muSpectre.gradient_integration functions. """ def setUp(self): self.lengths = np.array([2.4, 3.7, 4.1]) self.resolution = np.array([5, 3, 5]) self.norm_tol = 1e-8 def test_central_diff_derivative(self): """ Test of the central difference approximation by central_diff_derivative of the first derivative of a function on a grid. """ res = self.resolution * 15 lens = self.lengths for j in range(1, len(res)): d, dim, x_n, x_c, deriv, f_n, f_c, freqs = init_X_F_Chi(lens[:j], res[:j]) f_c = np.sin(2*np.pi/lens[:j] * x_c) for i in range(j): deriv[...,i,i] =2*np.pi/lens[i]*np.cos(2*np.pi* x_c[...,i]/lens[i]) approx_deriv = central_diff_derivative(f_c, d[:j+1], order=5) self.assertLess(np.linalg.norm(deriv-approx_deriv), self.norm_tol) def test_compute_wave_vectors(self): """ Test the construction of a wave vector grid by compute_wave_vectors for an arbitrary dimension. """ lens = [4, 6, 7] res = [3, 4, 5] Nx, Ny, Nz = res q_1d = lambda i: 2*np.pi/lens[i] * \ np.append(np.arange(res[i]-res[i]//2), -np.arange(1, res[i]//2 + 1)[::-1]) qx = q_1d(0) qy = q_1d(1) qz = q_1d(2) q = np.zeros(tuple(res) + (3,)) for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)): q[i,j,k,:] = np.array([qx[i], qy[j], qz[k]]) for n in range(1,4): comp_q = µ.gradient_integration.compute_wave_vectors(lens[:n], res[:n]) s = (np.s_[:],)*n + (0,)*(3-n) + (np.s_[:n],) self.assertLess(np.linalg.norm(comp_q - q[s]), self.norm_tol) def test_compute_grid(self): """ Test the function compute_grid which creates an orthogonal equally spaced grid of the given resolution and lengths. """ lens = self.lengths res = self.resolution d = np.array(lens)/np.array(res) grid_n = np.zeros(tuple(res+1) + (len(res),)) Nx, Ny, Nz = res+1 for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)): grid_n[i,j,k,:] = np.array([i*d[0], j*d[1], k*d[2]]) grid_c = (grid_n - d/2)[1:,1:,1:,:] for n in range(1,4): x_n, x_c = µ.gradient_integration.compute_grid(lens[:n], res[:n]) s = (np.s_[:],)*n + (0,)*(3-n) + (np.s_[:n],) self.assertLess(np.linalg.norm(x_c - grid_c[s]), self.norm_tol) self.assertLess(np.linalg.norm(x_n - grid_n[s]), self.norm_tol) def test_reshape_gradient(self): """ Test if reshape gradient transforms a flattend second order tensor in the right way to a shape resolution + [dim, dim]. """ lens = list(self.lengths) res = list(self.resolution) tol = 1e-5 formulation = µ.Formulation.finite_strain DelF = np.array([[0 , 0.01, 0.02], [0.03, 0 , 0.04], [0.05, 0.06, 0 ]]) one = np.eye(3,3) for n in range(2,4): sys = µ.Cell(res[:n], lens[:n], formulation) if n == 2: mat = µ.material.MaterialLinearElastic1_2d.make(sys, "material", 10, 0.3) if n == 3: mat = µ.material.MaterialLinearElastic1_3d.make(sys, "material", 10, 0.3) for pixel in sys: mat.add_pixel(pixel) solver = µ.solvers.SolverCG(sys, tol, maxiter=100, verbose=0) r = µ.solvers.newton_cg(sys, DelF[:n, :n], solver, tol, tol , verbose=0) grad = µ.gradient_integration.reshape_gradient(r.grad,list(res[:n])) grad_theo = (DelF[:n, :n] + one[:n, :n]).reshape((1,)*n+(n,n,)) self.assertEqual(grad.shape, tuple(res[:n])+(n,n,)) self.assertLess(np.linalg.norm(grad - grad_theo), self.norm_tol) def test_complement_periodically(self): """ Test the periodic reconstruction of an array. Lower left entries are added into the upper right part of the array. """ #1D grid scalars x_test = np.array([0,1,2,3]) x_test_p = np.array([0,1,2,3, 0]) x_p = µ.gradient_integration.complement_periodically(x_test, 1) self.assertLess(np.linalg.norm(x_p-x_test_p), self.norm_tol) #2D grid scalars x_test = np.array([[1,2,3,4], [5,6,7,8]]) x_test_p = np.array([[1,2,3,4,1], [5,6,7,8,5], [1,2,3,4,1]]) x_p = µ.gradient_integration.complement_periodically(x_test, 2) self.assertLess(np.linalg.norm(x_p-x_test_p), self.norm_tol) #2D grid vectors x_test = np.array([[[1,2,3] , [3,4,5]] , [[6,7,8] , [9,10,11]], [[12,13,14], [15,6,17]] ]) x_test_p = np.array([[[1,2,3] , [3,4,5] , [1,2,3]] , [[6,7,8] , [9,10,11], [6,7,8]] , [[12,13,14], [15,6,17], [12,13,14]], [[1,2,3] , [3,4,5] , [1,2,3]] ]) x_p = µ.gradient_integration.complement_periodically(x_test, 2) self.assertLess(np.linalg.norm(x_p-x_test_p), self.norm_tol) def test_get_integrator(self): """ Test if the right integrator is computed. """ #even grid lens_e = np.array([1,1,1]) res_e = np.array([2,2,2]) x_n_e, x_c_e = µ.gradient_integration.compute_grid(lens_e, res_e) freqs_e = µ.gradient_integration.compute_wave_vectors(lens_e, res_e) #odd grid lens_o = np.array([1,1]) res_o = np.array([3,3]) x_n_o, x_c_o = µ.gradient_integration.compute_grid(lens_o, res_o) delta_x = 1/3 freqs_o = µ.gradient_integration.compute_wave_vectors(lens_o, res_o) ### order=0 int_ana = 1j/(2*np.pi)*np.array([[[[ 0 , 0 , 0], [ 0 , 0 ,-1 ]] , [[ 0 ,-1 , 0], [ 0 ,-1/2,-1/2]]], [[[-1 , 0 , 0], [-1/2, 0 ,-1/2]] , [[-1/2,-1/2, 0], [-1/3,-1/3,-1/3]]] ]) dim,shape,integrator = µ.gradient_integration.\ get_integrator(x_c_e, freqs_e, order=0) self.assertEqual(dim, len(res_e)) self.assertEqual(shape, tuple(res_e)) self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol) ### order=1 #even grid int_ana = np.zeros(res_e.shape) dim,shape,integrator = µ.gradient_integration.\ get_integrator(x_c_e, freqs_e, order=1) self.assertEqual(dim, len(res_e)) self.assertEqual(shape, tuple(res_e)) self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol) #odd grid s = lambda q: np.sin(2*np.pi*q*delta_x) sq = lambda q: (np.sin(2*np.pi*np.array(q)*delta_x)**2).sum() int_ana = 1j * delta_x *\ np.array([[[ 0 , 0 ], [ 0 , s(1)/sq([0,1]) ], [ 0 , s(-1)/sq([0,-1]) ]], [[ s(1)/sq([1,0]) , 0 ], [ s(1)/sq([1,1]) , s(1)/sq([1,1]) ], [ s(1)/sq([1,-1]) , s(-1)/sq([1,-1]) ]], [[ s(-1)/sq([-1,0]) , 0 ], [ s(-1)/sq([-1,1]) , s(1)/sq([-1,1]) ], [ s(-1)/sq([-1,-1]) , s(-1)/sq([-1,-1]) ]]]) dim,shape,integrator = µ.gradient_integration.\ get_integrator(x_c_o, freqs_o, order=1) self.assertEqual(dim, len(res_o)) self.assertEqual(shape, tuple(res_o)) self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol) ### order=2 #even grid int_ana = np.zeros(res_e.shape) dim,shape,integrator = µ.gradient_integration.\ get_integrator(x_c_e, freqs_e, order=2) self.assertEqual(dim, len(res_e)) self.assertEqual(shape, tuple(res_e)) self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol) #odd grid lens_o = np.array([1,1]) res_o = np.array([3,3]) x_n, x_c = µ.gradient_integration.compute_grid(lens_o, res_o) delta_x = 1/3 freqs = µ.gradient_integration.compute_wave_vectors(lens_o, res_o) s = lambda q: 8*np.sin(2*np.pi*q*delta_x) + np.sin(2*2*np.pi*q*delta_x) sq = lambda q: ( (64*np.sin(2*np.pi*np.array(q)*delta_x)**2).sum() - (np.sin(2*2*np.pi*np.array(q)*delta_x)**2).sum() ) int_ana = 6 * 1j * delta_x *\ np.array([[[ 0 , 0 ], [ 0 , s(1)/sq([0,1]) ], [ 0 , s(-1)/sq([0,-1]) ]], [[ s(1)/sq([1,0]) , 0 ], [ s(1)/sq([1,1]) , s(1)/sq([1,1]) ], [ s(1)/sq([1,-1]) , s(-1)/sq([1,-1]) ]], [[ s(-1)/sq([-1,0]) , 0 ], [ s(-1)/sq([-1,1]) , s(1)/sq([-1,1]) ], [ s(-1)/sq([-1,-1]) , s(-1)/sq([-1,-1]) ]]]) dim,shape,integrator = µ.gradient_integration.\ get_integrator(x_c_o, freqs_o, order=2) self.assertEqual(dim, len(res_o)) self.assertEqual(shape, tuple(res_o)) self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol) def test_integrate_tensor_2(self): """ Test the correct integration of a second-rank tensor gradient field, like the deformation gradient. """ order = [1, 2] #list of higher order finite difference integration which #will be checked. ### cosinus, diagonal deformation gradient res = [15, 15, 14] lens = [7, 1.4, 3] d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res) for i in range(dim): F[:,:,:,i,i] = 0.8*np.pi/lens[i]*np.cos(4*np.pi* x_c[:,:,:,i]/lens[i]) Chi_n = 0.2 * np.sin(4*np.pi*x_n/lens) Chi_c = 0.2 * np.sin(4*np.pi*x_c/lens) # zeroth order correction placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs, staggered_grid=True, order=0) placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs, staggered_grid=False, order=0) self.assertLess(np.linalg.norm(Chi_n - placement_n), self.norm_tol) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) # higher order correction for n in order: tol_n = [1.334, 0.2299] #adjusted tolerances for node points F_c = central_diff_derivative(Chi_c, d, order=n) placement_n = µ.gradient_integration.integrate_tensor_2(F_c, x_n, freqs, staggered_grid=True, order=n) placement_c = µ.gradient_integration.integrate_tensor_2(F_c, x_c, freqs, staggered_grid=False,order=n) self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1]) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) ### cosinus, non-diagonal deformation gradient res = [15, 12, 11] lens = [8, 8, 8] d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res) F[:,:,:,0,0] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*x_c[:,:,:,0]) F[:,:,:,1,1] = 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*x_c[:,:,:,1]) F[:,:,:,2,2] = 2*np.pi/lens[2]*np.cos(2*np.pi/lens[2]*x_c[:,:,:,2]) F[:,:,:,1,0] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*x_c[:,:,:,0]) F[:,:,:,2,0] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*x_c[:,:,:,0]) for i in range(dim): Chi_c[:,:,:,i]= np.sin(2*np.pi*x_c[:,:,:,i]/lens[i]) \ + np.sin(2*np.pi*x_c[:,:,:,0]/lens[0]) Chi_n[:,:,:,i]= np.sin(2*np.pi*x_n[:,:,:,i]/lens[i]) \ + np.sin(2*np.pi*x_n[:,:,:,0]/lens[0]) # zeroth order correction placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs, staggered_grid=True, order=0) placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs, staggered_grid=False, order=0) self.assertLess(np.linalg.norm(Chi_n - placement_n), self.norm_tol) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) # higher order correction for n in order: tol_n = [2.563, 0.1544] #adjusted tolerances for node points F_c = central_diff_derivative(Chi_c, d, order=n) placement_n = µ.gradient_integration.integrate_tensor_2(F_c, x_n, freqs, staggered_grid=True, order=n) placement_c = µ.gradient_integration.integrate_tensor_2(F_c, x_c, freqs, staggered_grid=False,order=n) self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1]) self.assertLess(np.linalg.norm(Chi_c - placement_c), 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 = [13, 14, 11] lens = [4, 4, 4] d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res) for i in range(dim): F[:,:,:,i,i] = 32*x_c[:,:,:,i] -24*x_c[:,:,:,i]**2+4*x_c[:,:,:,i]**3 Chi_n = -128/15 + 16*x_n**2 -8*x_n**3 +x_n**4 Chi_c = -128/15 + 16*x_c**2 -8*x_c**3 +x_c**4 #subtract the mean of Chi_c, because the numeric integration is done to #give a zero mean fluctuating deformation field. mean_c = Chi_c.sum(axis=tuple(range(dim)))/ \ np.array(Chi_c.shape[:-1]).prod() Chi_n -= mean_c.reshape((1,)*dim + (dim,)) Chi_c -= mean_c.reshape((1,)*dim + (dim,)) # zeroth order correction placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs, staggered_grid=True, order=0) placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs, staggered_grid=False, order=0) self.assertLess(np.linalg.norm(Chi_n - placement_n), 0.19477) self.assertLess(np.linalg.norm(Chi_c - placement_c), 0.67355) # higher order correction for n in order: tol_n = [18.266, 2.9073] #adjusted tolerances for node points F_c = central_diff_derivative(Chi_c, d, order=n) placement_n = µ.gradient_integration.integrate_tensor_2(F_c, x_n, freqs, staggered_grid=True, order=n) placement_c = µ.gradient_integration.integrate_tensor_2(F_c, x_c, freqs, staggered_grid=False,order=n) self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1]) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) ### Realistic test: # shear of a two dimensional material with two different Young moduli. order_all = [0]+order #orders for which the test is run #initialize material structure res = [ 9, 21] #resolution lens = [ 9, 21] #lengths d, dim, x_n, x_c, _, _, _, freqs = init_X_F_Chi(lens, res) formulation = µ.Formulation.finite_strain Young = [10, 20] #Youngs modulus for each phase (soft, hard) Poisson = [0.3, 0.3] #Poissons ratio for each phase #geometry (two slabs stacked in y-direction with, #hight h (soft material) and hight res[1]-h (hard material)) h = res[1]//2 phase = np.zeros(tuple(res), dtype=int) phase[:, h:] = 1 phase = phase.flatten() cell = µ.Cell(res, lens, formulation) mat = µ.material.MaterialLinearElastic4_2d.make(cell, "material") for i, pixel in enumerate(cell): mat.add_pixel(pixel, Young[phase[i]], Poisson[phase[i]]) cell.initialise() DelF = np.array([[0 , 0.01], [0 , 0 ]]) # µSpectre solution solver = µ.solvers.SolverCG(cell, tol=1e-6, maxiter=100, verbose=0) result = µ.solvers.newton_cg(cell, DelF, solver, newton_tol=1e-6, equil_tol=1e-6, verbose=0) F = µ.gradient_integration.reshape_gradient(result.grad, res) fin_pos = {} #µSpectre computed center and node positions for all orders for n in order_all: placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs, staggered_grid=True, order=n) placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs, staggered_grid=False, order=n) fin_pos[str(n)+'_n'] = placement_n fin_pos[str(n)+'_c'] = placement_c # analytic solution, "placement_ana" (node and center) l_soft = d[1] * h #height soft material l_hard = d[1] * (res[1]-h) #height hard material Shear_modulus = np.array(Young) / (2 * (1+np.array(Poisson))) mean_shear_strain = 2*DelF[0,1] shear_strain_soft = (lens[1]*mean_shear_strain) / (l_soft + l_hard * Shear_modulus[0]/Shear_modulus[1]) shear_strain_hard = (lens[1]*mean_shear_strain) / (l_soft * Shear_modulus[1]/Shear_modulus[0] + l_hard) placement_ana_n = np.zeros(x_n.shape) placement_ana_c = np.zeros(x_c.shape) #x-coordinate #soft material placement_ana_n[:,:h+1,0] = shear_strain_soft/2 * x_n[:, :h+1, 1] placement_ana_c[:,:h ,0] = shear_strain_soft/2 * x_c[:, :h , 1] #hard material placement_ana_n[:,h+1:,0] =shear_strain_hard/2 * (x_n[:,h+1:,1]-l_soft)\ + shear_strain_soft/2 * l_soft placement_ana_c[:,h: ,0] =shear_strain_hard/2 * (x_c[:,h: ,1]-l_soft)\ + shear_strain_soft/2 * l_soft #y-coordinate placement_ana_n[:, :, 1] = 0 placement_ana_c[:, :, 1] = 0 #shift the analytic solution such that the average nonaffine deformation #is zero (integral of the nonaffine deformation gradient + N*const != 0) F_homo = (1./(np.prod(res)) * F.sum(axis=tuple(np.arange(dim))))\ .reshape((1,)*dim+(dim,)*2) #integration constant = integral of the nonaffine deformation gradient/N int_const = - ((placement_ana_c[:,:,0] - F_homo[:,:,0,1] * x_c[:,:,1]) .sum(axis=1))[0] / res[1] ana_sol_n = placement_ana_n + x_n + \ np.array([int_const, 0]).reshape((1,)*dim+(dim,)) ana_sol_c = placement_ana_c + x_c + \ np.array([int_const, 0]).reshape((1,)*dim+(dim,)) # check the numeric vs the analytic solution tol_n = [2.2112e-3, 1.3488e-3, 1.8124e-3] tol_c = [3.1095e-3, 3.2132e-2, 1.8989e-2] for n in order_all: norm_n = np.linalg.norm(fin_pos[str(n)+'_n'] - ana_sol_n) norm_c = np.linalg.norm(fin_pos[str(n)+'_c'] - ana_sol_c) self.assertLess(norm_n, tol_n[n]) self.assertLess(norm_c, tol_c[n]) def test_integrate_vector(self): """Test the integration of a first-rank tensor gradient field.""" order = [1,2] ### cosinus deformation gradient vector field res = [13, 14, 13] lens = [ 7, 4, 5] d, dim, x_n, x_c, df, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res, 1) for i in range(dim): df[:,:,:,i] = 0.8*np.pi/lens[i]*np.cos(4*np.pi*x_c[:,:,:,i]/lens[i]) Chi_n = 0.2 * np.sin(4*np.pi*x_n/lens).sum(axis=-1) Chi_c = 0.2 * np.sin(4*np.pi*x_c/lens).sum(axis=-1) # zeroth order correction placement_n = µ.gradient_integration.integrate_vector( df, x_n, freqs, staggered_grid=True, order=0) placement_c = µ.gradient_integration.integrate_vector( df, x_c, freqs, staggered_grid=False, order=0) self.assertLess(np.linalg.norm(Chi_n - placement_n), self.norm_tol) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) # higher order correction for n in order: tol_n = [1.404, 0.2882] #adjusted tolerances for node points df_c = central_diff_derivative(Chi_c, d, order=n, rank=0) placement_n = µ.gradient_integration.integrate_vector( df_c, x_n, freqs, staggered_grid=True, order=n) placement_c = µ.gradient_integration.integrate_vector( df_c, x_c, freqs, staggered_grid=False, order=n) self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1]) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) ### polynomial deformation gradient vector field # 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 = [12, 11, 13] lens = [4, 4, 4] d, dim, x_n, x_c, df, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res, 1) for i in range(dim): df[:,:,:,i] = 32*x_c[:,:,:,i]-24*x_c[:,:,:,i]**2+4*x_c[:,:,:,i]**3 Chi_n = (-128/15 + 16*x_n**2 -8*x_n**3 +x_n**4).sum(axis=-1) Chi_c = (-128/15 + 16*x_c**2 -8*x_c**3 +x_c**4).sum(axis=-1) #subtract the mean of Chi_c, because the numeric integration is done to #give a zero mean fluctuating deformation field. mean_c = Chi_c.sum() / np.array(Chi_c.shape).prod() Chi_n -= mean_c Chi_c -= mean_c # zeroth order correction placement_n = µ.gradient_integration.integrate_vector( df, x_n, freqs, staggered_grid=True, order=0) placement_c = µ.gradient_integration.integrate_vector( df, x_c, freqs, staggered_grid=False, order=0) self.assertLess(np.linalg.norm(Chi_n - placement_n), 0.20539) self.assertLess(np.linalg.norm(Chi_c - placement_c), 0.67380) # higher order correction for n in order: tol_n = [18.815, 3.14153] #adjusted tolerances for node points df_c = central_diff_derivative(Chi_c, d, order=n, rank=0) placement_n = µ.gradient_integration.integrate_vector( df_c, x_n, freqs, staggered_grid=True, order=n) placement_c = µ.gradient_integration.integrate_vector( df_c, x_c, freqs, staggered_grid=False, order=n) self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1]) self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol) def test_compute_placement(self): """Test the computation of placements and the original positions.""" ### shear of a homogeneous material ### res = [ 3, 11] #resolution lens = [10, 10] #lengths dim = len(res) #dimension x_n=µ.gradient_integration.compute_grid(np.array(lens),np.array(res))[0] ### finite strain formulation = µ.Formulation.finite_strain cell = µ.Cell(res, lens, formulation) mat = µ.material.MaterialLinearElastic1_2d.make(cell, "material", Young=10, Poisson=0.3) for pixel in cell: mat.add_pixel(pixel) cell.initialise() DelF = np.array([[0 , 0.05], [0 , 0 ]]) # analytic - placement_ana = np.zeros(x_n.shape) - placement_ana[:,:,0] = DelF[0,1]*x_n[:,:,1] + placement_ana = np.copy(x_n) + placement_ana[:,:,0] += DelF[0,1]*x_n[:,:,1] # µSpectre solution solver = µ.solvers.SolverCG(cell, tol=1e-6, maxiter=100, verbose=0) result = µ.solvers.newton_cg(cell, DelF, solver, newton_tol=1e-6, equil_tol=1e-6, verbose=0) for r in [result, result.grad]: #check input of result=OptimiseResult and result=np.ndarray placement, x = µ.gradient_integration.compute_placement( r, lens, res, order=0, formulation='finite_strain') self.assertLess(np.linalg.norm(placement_ana - placement), 1e-12) self.assertTrue((x_n == x).all())