223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

switched "Geom.materials" to "Geom.material"

This commit is contained in:
Philip Eisenlohr 2020-09-23 17:27:15 -04:00
parent b995f34834
commit 8c8db5b99f
7 changed files with 182 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
processing/pre/geom_grainGrowth.py
View File

@ -66,13 +66,13 @@ for name in filenames:
grid_original = geom.grid grid_original = geom.grid
damask.util.croak(geom) damask.util.croak(geom)
materials = np.tile(geom.materials,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1 material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1
grid = np.array(materials.shape) grid = np.array(material.shape)
# --- initialize support data --------------------------------------------------------------------- # --- initialize support data ---------------------------------------------------------------------
# store a copy of the initial material indices to find locations of immutable indices # store a copy of the initial material indices to find locations of immutable indices
materials_original = np.copy(materials) material_original = np.copy(material)
if not options.ndimage: if not options.ndimage:
X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]] X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]]
@ -88,14 +88,14 @@ for name in filenames:
for smoothIter in range(options.N): for smoothIter in range(options.N):
interfaceEnergy = np.zeros(materials.shape,dtype=np.float32) interfaceEnergy = np.zeros(material.shape,dtype=np.float32)
for i in (-1,0,1): for i in (-1,0,1):
for j in (-1,0,1): for j in (-1,0,1):
for k in (-1,0,1): for k in (-1,0,1):
# assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood) # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood)
interfaceEnergy = np.maximum(interfaceEnergy, interfaceEnergy = np.maximum(interfaceEnergy,
getInterfaceEnergy(materials,np.roll(np.roll(np.roll( getInterfaceEnergy(material,np.roll(np.roll(np.roll(
materials,i,axis=0), j,axis=1), k,axis=2))) 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 # 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, periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
@ -143,33 +143,33 @@ for name in filenames:
return_distances = False, return_distances = False,
return_indices = True) # want index of closest bulk grain return_indices = True) # want index of closest bulk grain
periodic_materials = np.tile(materials,(3,3,3))[grid[0]//2:-grid[0]//2, periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2, grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the geometry grid[2]//2:-grid[2]//2] # periodically extend the geometry
materials = periodic_materials[index[0], material = periodic_material[index[0],
index[1], index[1],
index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2, index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2, grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # extent grains into interface region grid[2]//2:-grid[2]//2] # extent grains into interface region
# replace immutable materials with closest mutable ones # replace immutable materials with closest mutable ones
index = ndimage.morphology.distance_transform_edt(np.in1d(materials,options.immutable).reshape(grid), index = ndimage.morphology.distance_transform_edt(np.in1d(material,options.immutable).reshape(grid),
return_distances = False, return_distances = False,
return_indices = True) return_indices = True)
materials = materials[index[0], material = material[index[0],
index[1], index[1],
index[2]] index[2]]
immutable = np.zeros(materials.shape, dtype=np.bool) immutable = np.zeros(material.shape, dtype=np.bool)
# find locations where immutable materials have been in original structure # find locations where immutable materials have been in original structure
for micro in options.immutable: for micro in options.immutable:
immutable += materials_original == micro immutable += material_original == micro
# undo any changes involving immutable materials # undo any changes involving immutable materials
materials = np.where(immutable, materials_original,materials) material = np.where(immutable, material_original,material)
damask.Geom(materials = materials[0:grid_original[0],0:grid_original[1],0:grid_original[2]], damask.Geom(material = material[0:grid_original[0],0:grid_original[1],0:grid_original[2]],
size = geom.size, size = geom.size,
origin = geom.origin, origin = geom.origin,
comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])], comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])],

12
processing/pre/mentat_spectralBox.py
View File

@ -100,7 +100,7 @@ def mesh(r,d):
#------------------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------------------
def material(): def materials():
return [\ return [\
"*new_mater standard", "*new_mater standard",
"*mater_option general:state:solid", "*mater_option general:state:solid",
@ -130,10 +130,10 @@ def geometry():
#------------------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------------------
def initial_conditions(materials): def initial_conditions(material):
elements = [] elements = []
element = 0 element = 0
for id in materials: for id in material:
element += 1 element += 1
if len(elements) < id: if len(elements) < id:
for i in range(id-len(elements)): for i in range(id-len(elements)):
@ -197,14 +197,14 @@ for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
materials = geom.materials.flatten(order='F') material = geom.material.flatten(order='F')
cmds = [\ cmds = [\
init(), init(),
mesh(geom.grid,geom.size), mesh(geom.grid,geom.size),
material(), materials(),
geometry(), geometry(),
initial_conditions(materials), initial_conditions(material),
'*identify_sets', '*identify_sets',
'*show_model', '*show_model',
'*redraw', '*redraw',

14
processing/pre/seeds_fromDistribution.py
View File

@ -101,8 +101,8 @@ class myThread (threading.Thread):
#--- evaluate current seeds file ------------------------------------------------------------------ #--- evaluate current seeds file ------------------------------------------------------------------
perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile) perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile)
myNmaterials = len(np.unique(perturbedGeom.materials)) myNmaterials = len(np.unique(perturbedGeom.material))
currentData=np.bincount(perturbedGeom.materials.ravel())[1:]/points currentData = np.bincount(perturbedGeom.material.ravel())[1:]/points
currentError=[] currentError=[]
currentHist=[] currentHist=[]
for i in range(nMaterials): # calculate the deviation in all bins per histogram for i in range(nMaterials): # calculate the deviation in all bins per histogram
@ -217,8 +217,8 @@ points = np.array(options.grid).prod().astype('float')
# ----------- calculate target distribution and bin edges # ----------- calculate target distribution and bin edges
targetGeom = damask.Geom.load_ASCII(os.path.splitext(os.path.basename(options.target))[0]+'.geom') targetGeom = damask.Geom.load_ASCII(os.path.splitext(os.path.basename(options.target))[0]+'.geom')
nMaterials = len(np.unique(targetGeom.materials)) nMaterials = len(np.unique(targetGeom.material))
targetVolFrac = np.bincount(targetGeom.materials.flatten())/targetGeom.grid.prod().astype(np.float) targetVolFrac = np.bincount(targetGeom.material.flatten())/targetGeom.grid.prod().astype(np.float)
target = [] target = []
for i in range(1,nMaterials+1): for i in range(1,nMaterials+1):
targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries
@ -244,10 +244,10 @@ initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
initialGeomVFile.seek(0) initialGeomVFile.seek(0)
initialGeom = damask.Geom.load_ASCII(initialGeomVFile) initialGeom = damask.Geom.load_ASCII(initialGeomVFile)
if len(np.unique(targetGeom.materials)) != nMaterials: if len(np.unique(targetGeom.material)) != nMaterials:
damask.util.croak('error. Material count mismatch') damask.util.croak('error. Material count mismatch')
initialData = np.bincount(initialGeom.materials.flatten())/points initialData = np.bincount(initialGeom.material.flatten())/points
for i in range(nMaterials): for i in range(nMaterials):
initialHist = np.histogram(initialData,bins=target[i]['bins'])[0] initialHist = np.histogram(initialData,bins=target[i]['bins'])[0]
target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum()) target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum())
@ -263,7 +263,7 @@ for i in range(nMaterials):
if options.maxseeds < 1: if options.maxseeds < 1:
maxSeeds = len(np.unique(initialGeom.materials)) maxSeeds = len(np.unique(initialGeom.material))
else: else:
maxSeeds = options.maxseeds maxSeeds = options.maxseeds

8
processing/pre/seeds_fromGeom.py
View File

@ -47,11 +47,11 @@ for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
materials = geom.materials.reshape((-1,1),order='F') material = geom.material.reshape((-1,1),order='F')
mask = np.logical_and(np.in1d(materials,options.whitelist,invert=False) if options.whitelist else \ mask = np.logical_and(np.in1d(material,options.whitelist,invert=False) if options.whitelist else \
np.full(geom.grid.prod(),True,dtype=bool), np.full(geom.grid.prod(),True,dtype=bool),
np.in1d(materials,options.blacklist,invert=True) if options.blacklist else \ np.in1d(material,options.blacklist,invert=True) if options.blacklist else \
np.full(geom.grid.prod(),True,dtype=bool)) np.full(geom.grid.prod(),True,dtype=bool))
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F') seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
@ -64,5 +64,5 @@ for name in filenames:
] ]
damask.Table(seeds[mask],{'pos':(3,)},comments)\ damask.Table(seeds[mask],{'pos':(3,)},comments)\
.add('material',materials[mask].astype(int))\ .add('material',material[mask].astype(int))\
.save(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds',legacy=True) .save(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds',legacy=True)

2
processing/pre/seeds_fromPokes.py
View File

@ -76,7 +76,7 @@ for name in filenames:
g[2] = k + offset[2] g[2] = k + offset[2]
g %= geom.grid g %= geom.grid
seeds[n,0:3] = (g+0.5)/geom.grid # normalize coordinates to box seeds[n,0:3] = (g+0.5)/geom.grid # normalize coordinates to box
seeds[n, 3] = geom.materials[g[0],g[1],g[2]] seeds[n, 3] = geom.material[g[0],g[1],g[2]]
if options.x: g[0] += 1 if options.x: g[0] += 1
if options.y: g[1] += 1 if options.y: g[1] += 1
n += 1 n += 1

234
python/damask/_geom.py
View File

@ -15,13 +15,13 @@ from . import grid_filters
class Geom: class Geom:
"""Geometry definition for grid solvers.""" """Geometry definition for grid solvers."""
def __init__(self,materials,size,origin=[0.0,0.0,0.0],comments=[]): def __init__(self,material,size,origin=[0.0,0.0,0.0],comments=[]):
""" """
New geometry definition from array of materials, size, and origin. New geometry definition from array of material, size, and origin.
Parameters Parameters
---------- ----------
materials : numpy.ndarray material : numpy.ndarray
Material index array (3D). Material index array (3D).
size : list or numpy.ndarray size : list or numpy.ndarray
Physical size of the geometry in meter. Physical size of the geometry in meter.
@ -31,16 +31,16 @@ class Geom:
Comment lines. Comment lines.
""" """
if len(materials.shape) != 3: if len(material.shape) != 3:
raise ValueError(f'Invalid materials shape {materials.shape}.') raise ValueError(f'Invalid material shape {material.shape}.')
elif materials.dtype not in np.sctypes['float'] + np.sctypes['int']: elif material.dtype not in np.sctypes['float'] + np.sctypes['int']:
raise TypeError(f'Invalid materials data type {materials.dtype}.') raise TypeError(f'Invalid material data type {material.dtype}.')
else: else:
self.materials = np.copy(materials) self.material = np.copy(material)
if self.materials.dtype in np.sctypes['float'] and \ if self.material.dtype in np.sctypes['float'] and \
np.all(self.materials == self.materials.astype(int).astype(float)): np.all(self.material == self.material.astype(int).astype(float)):
self.materials = self.materials.astype(int) self.material = self.material.astype(int)
if len(size) != 3 or any(np.array(size) <= 0): if len(size) != 3 or any(np.array(size) <= 0):
raise ValueError(f'Invalid size {size}.') raise ValueError(f'Invalid size {size}.')
@ -62,7 +62,7 @@ class Geom:
f'size x y z: {util.srepr(self.size, " x ")}', f'size x y z: {util.srepr(self.size, " x ")}',
f'origin x y z: {util.srepr(self.origin," ")}', f'origin x y z: {util.srepr(self.origin," ")}',
f'# materials: {self.N_materials}', f'# materials: {self.N_materials}',
f'max material: {np.nanmax(self.materials)}', f'max material: {np.nanmax(self.material)}',
]) ])
@ -103,21 +103,21 @@ class Geom:
message.append(util.delete(f'# materials: {other.N_materials}')) message.append(util.delete(f'# materials: {other.N_materials}'))
message.append(util.emph( f'# materials: { self.N_materials}')) message.append(util.emph( f'# materials: { self.N_materials}'))
if np.nanmax(other.materials) != np.nanmax(self.materials): if np.nanmax(other.material) != np.nanmax(self.material):
message.append(util.delete(f'max material: {np.nanmax(other.materials)}')) message.append(util.delete(f'max material: {np.nanmax(other.material)}'))
message.append(util.emph( f'max material: {np.nanmax( self.materials)}')) message.append(util.emph( f'max material: {np.nanmax( self.material)}'))
return util.return_message(message) return util.return_message(message)
@property @property
def grid(self): def grid(self):
return np.asarray(self.materials.shape) return np.asarray(self.material.shape)
@property @property
def N_materials(self): def N_materials(self):
return np.unique(self.materials).size return np.unique(self.material).size
@staticmethod @staticmethod
@ -160,7 +160,7 @@ class Geom:
else: else:
comments.append(line.strip()) comments.append(line.strip())
materials = np.empty(grid.prod()) # initialize as flat array material = np.empty(grid.prod()) # initialize as flat array
i = 0 i = 0
for line in content[header_length:]: for line in content[header_length:]:
items = line.split('#')[0].split() items = line.split('#')[0].split()
@ -172,16 +172,16 @@ class Geom:
abs(int(items[2])-int(items[0]))+1,dtype=float) abs(int(items[2])-int(items[0]))+1,dtype=float)
else: items = list(map(float,items)) else: items = list(map(float,items))
else: items = list(map(float,items)) else: items = list(map(float,items))
materials[i:i+len(items)] = items material[i:i+len(items)] = items
i += len(items) i += len(items)
if i != grid.prod(): if i != grid.prod():
raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}') raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}')
if not np.any(np.mod(materials,1) != 0.0): # no float present if not np.any(np.mod(material,1) != 0.0): # no float present
materials = materials.astype('int') material = material.astype('int')
return Geom(materials.reshape(grid,order='F'),size,origin,comments) return Geom(material.reshape(grid,order='F'),size,origin,comments)
@staticmethod @staticmethod
@ -200,9 +200,11 @@ class Geom:
comments = v.get_comments() comments = v.get_comments()
grid = np.array(v.vtk_data.GetDimensions())-1 grid = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
size = bbox[1] - bbox[0]
return Geom(v.get('material').reshape(grid,order='F'),size,bbox[0],comments=comments) return Geom(material = v.get('material').reshape(grid,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
@staticmethod @staticmethod
@ -243,17 +245,17 @@ class Geom:
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close() pool.close()
pool.join() pool.join()
materials = np.array(result.get()) material = np.array(result.get())
if periodic: if periodic:
materials = materials.reshape(grid*3) material = material.reshape(grid*3)
materials = materials[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] material = material[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else: else:
materials = materials.reshape(grid) material = material.reshape(grid)
return Geom(materials = materials+1, return Geom(material = material+1,
size = size, size = size,
comments = util.execution_stamp('Geom','from_Laguerre_tessellation'), comments = util.execution_stamp('Geom','from_Laguerre_tessellation'),
) )
@ -276,11 +278,11 @@ class Geom:
""" """
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,materials = KDTree.query(coords) devNull,material = KDTree.query(coords)
return Geom(materials = materials.reshape(grid)+1, return Geom(material = material.reshape(grid)+1,
size = size, size = size,
comments = util.execution_stamp('Geom','from_Voronoi_tessellation'), comments = util.execution_stamp('Geom','from_Voronoi_tessellation'),
) )
@ -311,10 +313,10 @@ class Geom:
plain = not compress plain = not compress
if plain: if plain:
format_string = '%g' if self.materials.dtype in np.sctypes['float'] else \ format_string = '%g' if self.material.dtype in np.sctypes['float'] else \
'%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.materials))))) '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material)))))
np.savetxt(fname, np.savetxt(fname,
self.materials.reshape([grid[0],np.prod(grid[1:])],order='F').T, self.material.reshape([grid[0],np.prod(grid[1:])],order='F').T,
header='\n'.join(header), fmt=format_string, comments='') header='\n'.join(header), fmt=format_string, comments='')
else: else:
try: try:
@ -325,7 +327,7 @@ class Geom:
compressType = None compressType = None
former = start = -1 former = start = -1
reps = 0 reps = 0
for current in self.materials.flatten('F'): for current in self.material.flatten('F'):
if abs(current - former) == 1 and (start - current) == reps*(former - current): if abs(current - former) == 1 and (start - current) == reps*(former - current):
compressType = 'to' compressType = 'to'
reps += 1 reps += 1
@ -370,7 +372,7 @@ class Geom:
""" """
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.add(self.materials.flatten(order='F'),'material') v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments) v.add_comments(self.comments)
v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress)
@ -402,7 +404,7 @@ class Geom:
0 gives octahedron (|x|^(2^0) + |y|^(2^0) + |z|^(2^0) < 1) 0 gives octahedron (|x|^(2^0) + |y|^(2^0) + |z|^(2^0) < 1)
1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1) 1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1)
fill : int, optional fill : int, optional
Fill value for primitive. Defaults to materials.max() + 1. Fill value for primitive. Defaults to material.max() + 1.
R : damask.Rotation, optional R : damask.Rotation, optional
Rotation of primitive. Defaults to no rotation. Rotation of primitive. Defaults to no rotation.
inverse : Boolean, optional inverse : Boolean, optional
@ -428,12 +430,12 @@ class Geom:
if periodic: # translate back to center if periodic: # translate back to center
mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2)) mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))
fill_ = np.full_like(self.materials,np.nanmax(self.materials)+1 if fill is None else fill) fill_ = np.full_like(self.material,np.nanmax(self.material)+1 if fill is None else fill)
return Geom(materials = np.where(np.logical_not(mask) if inverse else mask, self.materials,fill_), return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','add_primitive')], comments = self.comments+[util.execution_stamp('Geom','add_primitive')],
) )
@ -455,19 +457,19 @@ class Geom:
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.') raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
limits = [None,None] if reflect else [-2,0] limits = [None,None] if reflect else [-2,0]
ms = self.materials.copy() mat = self.material.copy()
if 'x' in directions: if 'x' in directions:
ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0) mat = np.concatenate([mat,mat[limits[0]:limits[1]:-1,:,:]],0)
if 'y' in directions: if 'y' in directions:
ms = np.concatenate([ms,ms[:,limits[0]:limits[1]:-1,:]],1) mat = np.concatenate([mat,mat[:,limits[0]:limits[1]:-1,:]],1)
if 'z' in directions: if 'z' in directions:
ms = np.concatenate([ms,ms[:,:,limits[0]:limits[1]:-1]],2) mat = np.concatenate([mat,mat[:,:,limits[0]:limits[1]:-1]],2)
return Geom(materials = ms, return Geom(material = mat,
size = self.size/self.grid*np.asarray(ms.shape), size = self.size/self.grid*np.asarray(mat.shape),
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','mirror')], comments = self.comments+[util.execution_stamp('Geom','mirror')],
) )
@ -486,12 +488,12 @@ class Geom:
if not set(directions).issubset(valid): if not set(directions).issubset(valid):
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.') raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
ms = np.flip(self.materials, (valid.index(d) for d in directions if d in valid)) mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
return Geom(materials = ms, return Geom(material = mat,
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','flip')], comments = self.comments+[util.execution_stamp('Geom','flip')],
) )
@ -507,17 +509,17 @@ class Geom:
Assume geometry to be periodic. Defaults to True. Assume geometry to be periodic. Defaults to True.
""" """
return Geom(materials = ndimage.interpolation.zoom( return Geom(material = ndimage.interpolation.zoom(
self.materials, self.material,
grid/self.grid, grid/self.grid,
output=self.materials.dtype, output=self.material.dtype,
order=0, order=0,
mode=('wrap' if periodic else 'nearest'), mode=('wrap' if periodic else 'nearest'),
prefilter=False prefilter=False
), ),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','scale')], comments = self.comments+[util.execution_stamp('Geom','scale')],
) )
@ -543,29 +545,29 @@ class Geom:
else: else:
return me return me
return Geom(materials = ndimage.filters.generic_filter( return Geom(material = ndimage.filters.generic_filter(
self.materials, self.material,
mostFrequent, mostFrequent,
size=(stencil if selection is None else stencil//2*2+1,)*3, size=(stencil if selection is None else stencil//2*2+1,)*3,
mode=('wrap' if periodic else 'nearest'), mode=('wrap' if periodic else 'nearest'),
extra_keywords=dict(selection=selection), extra_keywords=dict(selection=selection),
).astype(self.materials.dtype), ).astype(self.material.dtype),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','clean')], comments = self.comments+[util.execution_stamp('Geom','clean')],
) )
def renumber(self): def renumber(self):
"""Renumber sorted material indices to 1,...,N.""" """Renumber sorted material indices to 1,...,N."""
renumbered = np.empty(self.grid,dtype=self.materials.dtype) renumbered = np.empty(self.grid,dtype=self.material.dtype)
for i, oldID in enumerate(np.unique(self.materials)): for i, oldID in enumerate(np.unique(self.material)):
renumbered = np.where(self.materials == oldID, i+1, renumbered) renumbered = np.where(self.material == oldID, i+1, renumbered)
return Geom(materials = renumbered, return Geom(material = renumbered,
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','renumber')], comments = self.comments+[util.execution_stamp('Geom','renumber')],
) )
@ -578,32 +580,32 @@ class Geom:
R : damask.Rotation R : damask.Rotation
Rotation to apply to the geometry. Rotation to apply to the geometry.
fill : int or float, optional fill : int or float, optional
Material index to fill the corners. Defaults to materials.max() + 1. Material index to fill the corners. Defaults to material.max() + 1.
""" """
if fill is None: fill = np.nanmax(self.materials) + 1 if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.materials.dtype==np.float else int dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
Eulers = R.as_Eulers(degrees=True) Eulers = R.as_Eulers(degrees=True)
materials_in = self.materials.copy() material_in = self.material.copy()
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'') # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]): for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]):
materials_out = ndimage.rotate(materials_in,angle,axes,order=0, material_out = ndimage.rotate(material_in,angle,axes,order=0,
prefilter=False,output=dtype,cval=fill) prefilter=False,output=dtype,cval=fill)
if np.prod(materials_in.shape) == np.prod(materials_out.shape): if np.prod(material_in.shape) == np.prod(material_out.shape):
# avoid scipy interpolation errors for rotations close to multiples of 90° # avoid scipy interpolation errors for rotations close to multiples of 90°
materials_in = np.rot90(materials_in,k=np.rint(angle/90.).astype(int),axes=axes) material_in = np.rot90(material_in,k=np.rint(angle/90.).astype(int),axes=axes)
else: else:
materials_in = materials_out material_in = material_out
origin = self.origin-(np.asarray(materials_in.shape)-self.grid)*.5 * self.size/self.grid origin = self.origin-(np.asarray(material_in.shape)-self.grid)*.5 * self.size/self.grid
return Geom(materials = materials_in, return Geom(material = material_in,
size = self.size/self.grid*np.asarray(materials_in.shape), size = self.size/self.grid*np.asarray(material_in.shape),
origin = origin, origin = origin,
comments = self.comments+[util.execution_stamp('Geom','rotate')], comments = self.comments+[util.execution_stamp('Geom','rotate')],
) )
@ -618,12 +620,12 @@ class Geom:
offset : numpy.ndarray of shape (3) offset : numpy.ndarray of shape (3)
Offset (measured in grid points) from old to new geometry [0,0,0]. Offset (measured in grid points) from old to new geometry [0,0,0].
fill : int or float, optional fill : int or float, optional
Material index to fill the background. Defaults to materials.max() + 1. Material index to fill the background. Defaults to material.max() + 1.
""" """
if offset is None: offset = 0 if offset is None: offset = 0
if fill is None: fill = np.nanmax(self.materials) + 1 if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if int(fill) != fill or self.materials.dtype in np.sctypes['float'] else int dtype = float if int(fill) != fill or self.material.dtype in np.sctypes['float'] else int
canvas = np.full(self.grid if grid is None else grid,fill,dtype) canvas = np.full(self.grid if grid is None else grid,fill,dtype)
@ -632,35 +634,35 @@ class Geom:
ll = np.clip(-offset, 0,np.minimum( grid,self.grid-offset)) ll = np.clip(-offset, 0,np.minimum( grid,self.grid-offset))
ur = np.clip(-offset+self.grid,0,np.minimum( grid,self.grid-offset)) ur = np.clip(-offset+self.grid,0,np.minimum( grid,self.grid-offset))
canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.materials[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]] canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.material[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]]
return Geom(materials = canvas, return Geom(material = canvas,
size = self.size/self.grid*np.asarray(canvas.shape), size = self.size/self.grid*np.asarray(canvas.shape),
origin = self.origin+offset*self.size/self.grid, origin = self.origin+offset*self.size/self.grid,
comments = self.comments+[util.execution_stamp('Geom','canvas')], comments = self.comments+[util.execution_stamp('Geom','canvas')],
) )
def substitute(self,from_materials,to_materials): def substitute(self,from_material,to_material):
""" """
Substitute material indices. Substitute material indices.
Parameters Parameters
---------- ----------
from_materials : iterable of ints from_material : iterable of ints
Material indices to be substituted. Material indices to be substituted.
to_materials : iterable of ints to_material : iterable of ints
New material indices. New material indices.
""" """
substituted = self.materials.copy() substituted = self.material.copy()
for from_ms,to_ms in zip(from_materials,to_materials): for from_ms,to_ms in zip(from_material,to_material):
substituted[self.materials==from_ms] = to_ms substituted[self.material==from_ms] = to_ms
return Geom(materials = substituted, return Geom(material = substituted,
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','substitute')], comments = self.comments+[util.execution_stamp('Geom','substitute')],
) )
@ -679,7 +681,7 @@ class Geom:
Defaults to 1. Defaults to 1.
offset : int, optional offset : int, optional
Offset (positive or negative) to tag material indices, Offset (positive or negative) to tag material indices,
defaults to materials.max() + 1. defaults to material.max() + 1.
trigger : list of ints, optional trigger : list of ints, optional
List of material indices that trigger a change. List of material indices that trigger a change.
Defaults to [], meaning that any different neighbor triggers a change. Defaults to [], meaning that any different neighbor triggers a change.
@ -698,15 +700,15 @@ class Geom:
trigger = list(trigger) trigger = list(trigger)
return np.any(np.in1d(stencil,np.array(trigger))) return np.any(np.in1d(stencil,np.array(trigger)))
offset_ = np.nanmax(self.materials) if offset is None else offset offset_ = np.nanmax(self.material) if offset is None else offset
mask = ndimage.filters.generic_filter(self.materials, mask = ndimage.filters.generic_filter(self.material,
tainted_neighborhood, tainted_neighborhood,
size=1+2*vicinity, size=1+2*vicinity,
mode='wrap' if periodic else 'nearest', mode='wrap' if periodic else 'nearest',
extra_keywords={'trigger':trigger}) extra_keywords={'trigger':trigger})
return Geom(materials = np.where(mask, self.materials + offset_,self.materials), return Geom(material = np.where(mask, self.material + offset_,self.material),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')], comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')],
) )

52
python/tests/test_Geom.py
View File

@ -11,8 +11,8 @@ from damask import util
def geom_equal(a,b): def geom_equal(a,b):
return np.all(a.materials == b.materials) and \ return np.all(a.material == b.material) and \
np.all(a.grid == b.grid) and \ np.all(a.grid == b.grid) and \
np.allclose(a.size, b.size) and \ np.allclose(a.size, b.size) and \
str(a.diff(b)) == str(b.diff(a)) str(a.diff(b)) == str(b.diff(a))
@ -38,7 +38,7 @@ class TestGeom:
def test_diff_not_equal(self,default): def test_diff_not_equal(self,default):
new = Geom(default.materials[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified']) new = Geom(default.material[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified'])
assert str(default.diff(new)) != '' assert str(default.diff(new)) != ''
@ -93,28 +93,28 @@ class TestGeom:
def test_invalid_size(self,default): def test_invalid_size(self,default):
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom(default.materials[1:,1:,1:], Geom(default.material[1:,1:,1:],
size=np.ones(2)) size=np.ones(2))
def test_invalid_origin(self,default): def test_invalid_origin(self,default):
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom(default.materials[1:,1:,1:], Geom(default.material[1:,1:,1:],
size=np.ones(3), size=np.ones(3),
origin=np.ones(4)) origin=np.ones(4))
def test_invalid_materials_shape(self,default): def test_invalid_materials_shape(self,default):
materials = np.ones((3,3)) material = np.ones((3,3))
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom(materials, Geom(material,
size=np.ones(3)) size=np.ones(3))
def test_invalid_materials_type(self,default): def test_invalid_materials_type(self,default):
materials = np.random.randint(1,300,(3,4,5))==1 material = np.random.randint(1,300,(3,4,5))==1
with pytest.raises(TypeError): with pytest.raises(TypeError):
Geom(materials) Geom(material)
@pytest.mark.parametrize('directions,reflect',[ @pytest.mark.parametrize('directions,reflect',[
@ -205,10 +205,10 @@ class TestGeom:
def test_renumber(self,default): def test_renumber(self,default):
materials = default.materials.copy() material = default.material.copy()
for m in np.unique(materials): for m in np.unique(material):
materials[materials==m] = materials.max() + np.random.randint(1,30) material[material==m] = material.max() + np.random.randint(1,30)
modified = Geom(materials, modified = Geom(material,
default.size, default.size,
default.origin) default.origin)
assert not geom_equal(modified,default) assert not geom_equal(modified,default)
@ -218,13 +218,13 @@ class TestGeom:
def test_substitute(self,default): def test_substitute(self,default):
offset = np.random.randint(1,500) offset = np.random.randint(1,500)
modified = Geom(default.materials + offset, modified = Geom(default.material + offset,
default.size, default.size,
default.origin) default.origin)
assert not geom_equal(modified,default) assert not geom_equal(modified,default)
assert geom_equal(default, assert geom_equal(default,
modified.substitute(np.arange(default.materials.max())+1+offset, modified.substitute(np.arange(default.material.max())+1+offset,
np.arange(default.materials.max())+1)) np.arange(default.material.max())+1))
@pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]), @pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
@ -251,7 +251,7 @@ class TestGeom:
grid = default.grid grid = default.grid
grid_add = np.random.randint(0,30,(3)) grid_add = np.random.randint(0,30,(3))
modified = default.canvas(grid + grid_add) modified = default.canvas(grid + grid_add)
assert np.all(modified.materials[:grid[0],:grid[1],:grid[2]] == default.materials) assert np.all(modified.material[:grid[0],:grid[1],:grid[2]] == default.material)
@pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random()*8), @pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random()*8),
@ -271,7 +271,7 @@ class TestGeom:
s = np.random.random(3)+.5 s = np.random.random(3)+.5
G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,exponent) G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,exponent)
G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,exponent) G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,exponent)
assert np.count_nonzero(G_1.materials!=2) == np.count_nonzero(G_2.materials!=2) assert np.count_nonzero(G_1.material!=2) == np.count_nonzero(G_2.material!=2)
@pytest.mark.parametrize('center',[np.random.randint(4,10,(3)), @pytest.mark.parametrize('center',[np.random.randint(4,10,(3)),
@ -308,14 +308,14 @@ class TestGeom:
geom = Geom(m,np.random.rand(3)).vicinity_offset(vicinity,offset,trigger=trigger) geom = Geom(m,np.random.rand(3)).vicinity_offset(vicinity,offset,trigger=trigger)
assert np.all(m2==geom.materials) assert np.all(m2==geom.material)
@pytest.mark.parametrize('periodic',[True,False]) @pytest.mark.parametrize('periodic',[True,False])
def test_vicinity_offset_invariant(self,default,periodic): def test_vicinity_offset_invariant(self,default,periodic):
offset = default.vicinity_offset(trigger=[default.materials.max()+1, offset = default.vicinity_offset(trigger=[default.material.max()+1,
default.materials.min()-1]) default.material.min()-1])
assert np.all(offset.materials==default.materials) assert np.all(offset.material==default.material)
@pytest.mark.parametrize('periodic',[True,False]) @pytest.mark.parametrize('periodic',[True,False])
@ -338,7 +338,7 @@ class TestGeom:
ms = np.random.randint(1, N_seeds+1) ms = np.random.randint(1, N_seeds+1)
weights[ms-1] = np.random.random() weights[ms-1] = np.random.random()
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5) Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5)
assert np.all(Laguerre.materials == ms) assert np.all(Laguerre.material == ms)
@pytest.mark.parametrize('approach',['Laguerre','Voronoi']) @pytest.mark.parametrize('approach',['Laguerre','Voronoi'])
@ -346,10 +346,10 @@ class TestGeom:
grid = np.random.randint(5,10,3)*2 grid = np.random.randint(5,10,3)*2
size = grid.astype(np.float) size = grid.astype(np.float)
seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5]))) seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5])))
materials = np.ones(grid) material = np.ones(grid)
materials[:,grid[1]//2:,:] = 2 material[:,grid[1]//2:,:] = 2
if approach == 'Laguerre': if approach == 'Laguerre':
geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5) geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5)
elif approach == 'Voronoi': elif approach == 'Voronoi':
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5) geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5)
assert np.all(geom.materials == materials) assert np.all(geom.material == material)