diff --git a/processing/legacy/geom_grainGrowth.py b/processing/legacy/geom_grainGrowth.py deleted file mode 100755 index d969b415c..000000000 --- a/processing/legacy/geom_grainGrowth.py +++ /dev/null @@ -1,178 +0,0 @@ -#!/usr/bin/env python3 - -import os -import sys -from io import StringIO -from optparse import OptionParser - -import numpy as np -from scipy import ndimage - -import damask - - -scriptName = os.path.splitext(os.path.basename(__file__))[0] -scriptID = ' '.join([scriptName,damask.version]) - - -getInterfaceEnergy = lambda A,B: np.float32((A != B)*1.0) # 1.0 if A & B are distinct, 0.0 otherwise -struc = ndimage.generate_binary_structure(3,1) # 3D von Neumann neighborhood - - -#-------------------------------------------------------------------------------------------------- -# MAIN -#-------------------------------------------------------------------------------------------------- - -parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [geomfile(s)]', description = """ -Smoothen interface roughness by simulated curvature flow. -This is achieved by the diffusion of each initially sharply bounded grain volume within the periodic domain -up to a given distance 'd' voxels. -The final geometry is assembled by selecting at each voxel that grain index for which the concentration remains largest. -References 10.1073/pnas.1111557108 (10.1006/jcph.1994.1105) - -""", version = scriptID) - -parser.add_option('-d', '--distance', - dest = 'd', - type = 'float', metavar = 'float', - help = 'diffusion distance in voxels [%default]') -parser.add_option('-N', '--iterations', - dest = 'N', - type = 'int', metavar = 'int', - help = 'curvature flow iterations [%default]') -parser.add_option('-i', '--immutable', - action = 'extend', dest = 'immutable', metavar = '', - help = 'list of immutable material indices') -parser.add_option('--ndimage', - dest = 'ndimage', action='store_true', - help = 'use ndimage.gaussian_filter in lieu of explicit FFT') - -parser.set_defaults(d = 1, - N = 1, - immutable = [], - ndimage = False, - ) - -(options, filenames) = parser.parse_args() - -options.immutable = list(map(int,options.immutable)) - - -if filenames == []: filenames = [None] - -for name in filenames: - damask.util.report(scriptName,name) - - geom = damask.Grid.load(StringIO(''.join(sys.stdin.read())) if name is None else name) - - grid_original = geom.cells - damask.util.croak(geom) - material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1 - grid = np.array(material.shape) - -# --- initialize support data --------------------------------------------------------------------- - -# store a copy of the initial material indices to find locations of immutable indices - material_original = np.copy(material) - - if not options.ndimage: - X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]] - - # Calculates gaussian weights for simulating 3d diffusion - gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d),dtype=np.float32) \ - /np.power(2.0*np.pi*options.d*options.d,(3.0 - np.count_nonzero(grid_original == 1))/2.,dtype=np.float32) - - gauss[:,:,:grid[2]//2:-1] = gauss[:,:,1:(grid[2]+1)//2] # trying to cope with uneven (odd) grid size - gauss[:,:grid[1]//2:-1,:] = gauss[:,1:(grid[1]+1)//2,:] - gauss[:grid[0]//2:-1,:,:] = gauss[1:(grid[0]+1)//2,:,:] - gauss = np.fft.rfftn(gauss).astype(np.complex64) - - for smoothIter in range(options.N): - - interfaceEnergy = np.zeros(material.shape,dtype=np.float32) - for i in (-1,0,1): - for j in (-1,0,1): - for k in (-1,0,1): - # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood) - interfaceEnergy = np.maximum(interfaceEnergy, - getInterfaceEnergy(material,np.roll(np.roll(np.roll( - material,i,axis=0), j,axis=1), k,axis=2))) - - # periodically extend interfacial energy array by half a grid size in positive and negative directions - periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] - - # transform bulk volume (i.e. where interfacial energy remained zero), store index of closest boundary voxel - index = ndimage.morphology.distance_transform_edt(periodic_interfaceEnergy == 0., - return_distances = False, - return_indices = True) - - # want array index of nearest voxel on periodically extended boundary - periodic_bulkEnergy = periodic_interfaceEnergy[index[0], - index[1], - index[2]].reshape(2*grid) # fill bulk with energy of nearest interface - - if options.ndimage: - periodic_diffusedEnergy = ndimage.gaussian_filter( - np.where(ndimage.morphology.binary_dilation(periodic_interfaceEnergy > 0., - structure = struc, - iterations = int(round(options.d*2.))-1, # fat boundary - ), - periodic_bulkEnergy, # ...and zero everywhere else - 0.), - sigma = options.d) - else: - diffusedEnergy = np.fft.irfftn(np.fft.rfftn( - np.where( - ndimage.morphology.binary_dilation(interfaceEnergy > 0., - structure = struc, - iterations = int(round(options.d*2.))-1),# fat boundary - periodic_bulkEnergy[grid[0]//2:-grid[0]//2, # retain filled energy on fat boundary... - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2], # ...and zero everywhere else - 0.)).astype(np.complex64) * - gauss).astype(np.float32) - - periodic_diffusedEnergy = np.tile(diffusedEnergy,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # periodically extend the smoothed bulk energy - - - # transform voxels close to interface region - index = ndimage.morphology.distance_transform_edt(periodic_diffusedEnergy >= 0.95*np.amax(periodic_diffusedEnergy), - return_distances = False, - return_indices = True) # want index of closest bulk grain - - periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # periodically extend the geometry - - material = periodic_material[index[0], - index[1], - index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # extent grains into interface region - - # replace immutable materials with closest mutable ones - index = ndimage.morphology.distance_transform_edt(np.in1d(material,options.immutable).reshape(grid), - return_distances = False, - return_indices = True) - material = material[index[0], - index[1], - index[2]] - - immutable = np.zeros(material.shape, dtype=np.bool) - # find locations where immutable materials have been in original structure - for micro in options.immutable: - immutable += material_original == micro - - # undo any changes involving immutable materials - material = np.where(immutable, material_original,material) - - damask.Grid(material = material[0:grid_original[0],0:grid_original[1],0:grid_original[2]], - size = geom.size, - origin = geom.origin, - comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])], - )\ - .save(sys.stdout if name is None else name) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 2a32adea4..bb4a14ed7 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -287,7 +287,7 @@ class TestOrientation: @pytest.mark.parametrize('family',crystal_families) @pytest.mark.parametrize('proper',[True,False]) def test_in_SST(self,family,proper): - assert Orientation(family=family).in_SST(np.zeros(3),proper) + assert Orientation(family=family).in_SST(np.zeros(3),proper) @pytest.mark.parametrize('function',['in_SST','IPF_color']) def test_invalid_argument(self,function): diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 1da56dedf..71a4ce802 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -367,13 +367,13 @@ class TestResult: @pytest.mark.parametrize('mode',['cell','node']) def test_coordinates(self,default,mode): - if mode == 'cell': - a = grid_filters.coordinates0_point(default.cells,default.size,default.origin) - b = default.coordinates0_point.reshape(tuple(default.cells)+(3,),order='F') - elif mode == 'node': - a = grid_filters.coordinates0_node(default.cells,default.size,default.origin) - b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F') - assert np.allclose(a,b) + if mode == 'cell': + a = grid_filters.coordinates0_point(default.cells,default.size,default.origin) + b = default.coordinates0_point.reshape(tuple(default.cells)+(3,),order='F') + elif mode == 'node': + a = grid_filters.coordinates0_node(default.cells,default.size,default.origin) + b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F') + assert np.allclose(a,b) @pytest.mark.parametrize('output',['F','*',['P'],['P','F']],ids=range(4)) @pytest.mark.parametrize('fname',['12grains6x7x8_tensionY.hdf5'],ids=range(1)) @@ -421,7 +421,7 @@ class TestResult: def test_XDMF_datatypes(self,tmp_path,single_phase,update,ref_path): for shape in [('scalar',()),('vector',(3,)),('tensor',(3,3)),('matrix',(12,))]: for dtype in ['f4','f8','i1','i2','i4','i8','u1','u2','u4','u8']: - single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}') + single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}') fname = os.path.splitext(os.path.basename(single_phase.fname))[0]+'.xdmf' os.chdir(tmp_path) single_phase.export_XDMF() diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a431bc64b..fda986c2f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -1076,19 +1076,19 @@ class TestRotation: def test_from_fiber_component(self,N,sigma): p = [] for run in range(5): - alpha = np.random.random()*2*np.pi,np.arccos(np.random.random()) - beta = np.random.random()*2*np.pi,np.arccos(np.random.random()) + alpha = np.random.random()*2*np.pi,np.arccos(np.random.random()) + beta = np.random.random()*2*np.pi,np.arccos(np.random.random()) - f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])]) - f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )]) - ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S))) - n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system + f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])]) + f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )]) + ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S))) + n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system - o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False) - angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1)) - dist = np.array(angles) * (np.random.randint(0,2,N)*2-1) + o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False) + angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1)) + dist = np.array(angles) * (np.random.randint(0,2,N)*2-1) - p.append(stats.normaltest(dist)[1]) + p.append(stats.normaltest(dist)[1]) sigma_out = np.degrees(np.std(dist)) p = np.average(p) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index e59409a20..9f0dcc7cf 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -173,11 +173,11 @@ class TestVTK: polyData = VTK.from_poly_data(points) polyData.add(points,'coordinates') if update: - polyData.save(ref_path/'polyData') + polyData.save(ref_path/'polyData') else: - reference = VTK.load(ref_path/'polyData.vtp') - assert polyData.__repr__() == reference.__repr__() and \ - np.allclose(polyData.get('coordinates'),points) + reference = VTK.load(ref_path/'polyData.vtp') + assert polyData.__repr__() == reference.__repr__() and \ + np.allclose(polyData.get('coordinates'),points) @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA') def test_compare_reference_rectilinearGrid(self,update,ref_path,tmp_path): @@ -189,8 +189,8 @@ class TestVTK: rectilinearGrid.add(np.ascontiguousarray(c),'cell') rectilinearGrid.add(np.ascontiguousarray(n),'node') if update: - rectilinearGrid.save(ref_path/'rectilinearGrid') + rectilinearGrid.save(ref_path/'rectilinearGrid') else: - reference = VTK.load(ref_path/'rectilinearGrid.vtr') - assert rectilinearGrid.__repr__() == reference.__repr__() and \ - np.allclose(rectilinearGrid.get('cell'),c) + reference = VTK.load(ref_path/'rectilinearGrid.vtr') + assert rectilinearGrid.__repr__() == reference.__repr__() and \ + np.allclose(rectilinearGrid.get('cell'),c)