DAMASK_EICMD/python/damask/_geom.py

721 lines
28 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

import copy
import multiprocessing as mp
from functools import partial
import numpy as np
from scipy import ndimage,spatial
from . import environment
from . import Rotation
from . import VTK
from . import util
from . import grid_filters
class Geom:
"""Geometry definition for grid solvers."""
def __init__(self,material,size,origin=[0.0,0.0,0.0],comments=[]):
"""
New geometry definition from array of material, size, and origin.
Parameters
----------
material : numpy.ndarray
Material index array (3D).
size : list or numpy.ndarray
Physical size of the geometry in meter.
origin : list or numpy.ndarray, optional
Physical origin of the geometry in meter.
comments : list of str, optional
Comment lines.
"""
if len(material.shape) != 3:
raise ValueError(f'Invalid material shape {material.shape}.')
elif material.dtype not in np.sctypes['float'] + np.sctypes['int']:
raise TypeError(f'Invalid material data type {material.dtype}.')
else:
self.material = np.copy(material)
if self.material.dtype in np.sctypes['float'] and \
np.all(self.material == self.material.astype(int).astype(float)):
self.material = self.material.astype(int)
if len(size) != 3 or any(np.array(size) <= 0):
raise ValueError(f'Invalid size {size}.')
else:
self.size = np.array(size)
if len(origin) != 3:
raise ValueError(f'Invalid origin {origin}.')
else:
self.origin = np.array(origin)
self.comments = [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
def __repr__(self):
"""Basic information on geometry definition."""
return util.srepr([
f'grid a b c: {util.srepr(self.grid, " x ")}',
f'size x y z: {util.srepr(self.size, " x ")}',
f'origin x y z: {util.srepr(self.origin," ")}',
f'# materials: {self.N_materials}',
f'max material: {np.nanmax(self.material)}',
])
def __copy__(self):
"""Copy geometry."""
return copy.deepcopy(self)
def copy(self):
"""Copy geometry."""
return self.__copy__()
def diff(self,other):
"""
Report property differences of self relative to other.
Parameters
----------
other : Geom
Geometry to compare self against.
"""
message = []
if np.any(other.grid != self.grid):
message.append(util.delete(f'grid a b c: {util.srepr(other.grid," x ")}'))
message.append(util.emph( f'grid a b c: {util.srepr( self.grid," x ")}'))
if not np.allclose(other.size,self.size):
message.append(util.delete(f'size x y z: {util.srepr(other.size," x ")}'))
message.append(util.emph( f'size x y z: {util.srepr( self.size," x ")}'))
if not np.allclose(other.origin,self.origin):
message.append(util.delete(f'origin x y z: {util.srepr(other.origin," ")}'))
message.append(util.emph( f'origin x y z: {util.srepr( self.origin," ")}'))
if other.N_materials != self.N_materials:
message.append(util.delete(f'# materials: {other.N_materials}'))
message.append(util.emph( f'# materials: { self.N_materials}'))
if np.nanmax(other.material) != np.nanmax(self.material):
message.append(util.delete(f'max material: {np.nanmax(other.material)}'))
message.append(util.emph( f'max material: {np.nanmax( self.material)}'))
return util.return_message(message)
@property
def grid(self):
return np.asarray(self.material.shape)
@property
def N_materials(self):
return np.unique(self.material).size
@staticmethod
def load_ASCII(fname):
"""
Read a geom file.
Parameters
----------
fname : str or file handle
Geometry file to read.
"""
try:
f = open(fname)
except TypeError:
f = fname
f.seek(0)
try:
header_length,keyword = f.readline().split()[:2]
header_length = int(header_length)
except ValueError:
header_length,keyword = (-1, 'invalid')
if not keyword.startswith('head') or header_length < 3:
raise TypeError('Header length information missing or invalid')
content = f.readlines()
comments = []
for i,line in enumerate(content[:header_length]):
items = line.split('#')[0].lower().strip().split()
key = items[0] if items else ''
if key == 'grid':
grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']])
elif key == 'size':
size = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
elif key == 'origin':
origin = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
else:
comments.append(line.strip())
material = np.empty(grid.prod()) # initialize as flat array
i = 0
for line in content[header_length:]:
items = line.split('#')[0].split()
if len(items) == 3:
if items[1].lower() == 'of':
items = np.ones(int(items[0]))*float(items[2])
elif items[1].lower() == 'to':
items = np.linspace(int(items[0]),int(items[2]),
abs(int(items[2])-int(items[0]))+1,dtype=float)
else: items = list(map(float,items))
else: items = list(map(float,items))
material[i:i+len(items)] = items
i += len(items)
if i != grid.prod():
raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}')
if not np.any(np.mod(material,1) != 0.0): # no float present
material = material.astype('int')
return Geom(material.reshape(grid,order='F'),size,origin,comments)
@staticmethod
def load(fname):
"""
Read a VTK rectilinear grid.
Parameters
----------
fname : str or or pathlib.Path
Geometry file to read.
Valid extension is .vtr, it will be appended if not given.
"""
v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
comments = v.get_comments()
grid = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
return Geom(material = v.get('material').reshape(grid,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
@staticmethod
def _find_closest_seed(seeds, weights, point):
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@staticmethod
def from_Laguerre_tessellation(grid,size,seeds,weights,material=None,periodic=True):
"""
Generate geometry from Laguerre tessellation.
Parameters
----------
grid : int numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
Physical size of the geometry in meter.
seeds : numpy.ndarray of shape (:,3)
Position of the seed points in meter. All points need to lay within the box.
weights : numpy.ndarray of shape (seeds.shape[0])
Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered.
periodic : Boolean, optional
Perform a periodic tessellation. Defaults to True.
"""
if periodic:
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3)
else:
weights_p = weights
seeds_p = seeds
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
pool = mp.Pool(processes = int(environment.options['DAMASK_NUM_THREADS']))
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()
material_ = np.array(result.get())
if periodic:
material_ = material_.reshape(grid*3)
material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else:
material_ = material_.reshape(grid)
return Geom(material = material_+1 if material is None else material[material_],
size = size,
comments = util.execution_stamp('Geom','from_Laguerre_tessellation'),
)
@staticmethod
def from_Voronoi_tessellation(grid,size,seeds,material=None,periodic=True):
"""
Generate geometry from Voronoi tessellation.
Parameters
----------
grid : int numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
Physical size of the geometry in meter.
seeds : numpy.ndarray of shape (:,3)
Position of the seed points in meter. All points need to lay within the box.
material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered.
periodic : Boolean, optional
Perform a periodic tessellation. Defaults to True.
"""
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,material_ = KDTree.query(coords)
return Geom(material = (material_+1 if material is None else material[material_]).reshape(grid),
size = size,
comments = util.execution_stamp('Geom','from_Voronoi_tessellation'),
)
def save_ASCII(self,fname,compress=None):
"""
Writes a geom file.
Parameters
----------
fname : str or file handle
Geometry file to write with extension '.geom'.
compress : bool, optional
Compress geometry with 'x of y' and 'a to b'.
"""
header = [f'{len(self.comments)+4} header'] + self.comments \
+ ['grid a {} b {} c {}'.format(*self.grid),
'size x {} y {} z {}'.format(*self.size),
'origin x {} y {} z {}'.format(*self.origin),
'homogenization 1',
]
grid = self.grid
if compress is None:
plain = grid.prod()/self.N_materials < 250
else:
plain = not compress
if plain:
format_string = '%g' if self.material.dtype in np.sctypes['float'] else \
'%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material)))))
np.savetxt(fname,
self.material.reshape([grid[0],np.prod(grid[1:])],order='F').T,
header='\n'.join(header), fmt=format_string, comments='')
else:
try:
f = open(fname,'w')
except TypeError:
f = fname
compressType = None
former = start = -1
reps = 0
for current in self.material.flatten('F'):
if abs(current - former) == 1 and (start - current) == reps*(former - current):
compressType = 'to'
reps += 1
elif current == former and start == former:
compressType = 'of'
reps += 1
else:
if compressType is None:
f.write('\n'.join(header)+'\n')
elif compressType == '.':
f.write(f'{former}\n')
elif compressType == 'to':
f.write(f'{start} to {former}\n')
elif compressType == 'of':
f.write(f'{reps} of {former}\n')
compressType = '.'
start = current
reps = 1
former = current
if compressType == '.':
f.write(f'{former}\n')
elif compressType == 'to':
f.write(f'{start} to {former}\n')
elif compressType == 'of':
f.write(f'{reps} of {former}\n')
def save(self,fname,compress=True):
"""
Generates vtk rectilinear grid.
Parameters
----------
fname : str, optional
Filename to write. If no file is given, a string is returned.
Valid extension is .vtr, it will be appended if not given.
compress : bool, optional
Compress with zlib algorithm. Defaults to True.
"""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments)
v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress)
def show(self):
"""Show on screen."""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.show()
def add_primitive(self,dimension,center,exponent,
fill=None,R=Rotation(),inverse=False,periodic=True):
"""
Inserts a primitive geometric object at a given position.
Parameters
----------
dimension : int or float numpy.ndarray of shape(3)
Dimension (diameter/side length) of the primitive. If given as
integers, grid point locations (cell centers) are addressed.
If given as floats, coordinates are addressed.
center : int or float numpy.ndarray of shape(3)
Center of the primitive. If given as integers, grid point
locations (cell centers) are addressed.
If given as floats, coordinates are addressed.
exponent : numpy.ndarray of shape(3) or float
Exponents for the three axes.
0 gives octahedron (ǀxǀ^(2^0) + ǀyǀ^(2^0) + ǀzǀ^(2^0) < 1)
1 gives sphere (ǀxǀ^(2^1) + ǀyǀ^(2^1) + ǀzǀ^(2^1) < 1)
fill : int, optional
Fill value for primitive. Defaults to material.max() + 1.
R : damask.Rotation, optional
Rotation of primitive. Defaults to no rotation.
inverse : Boolean, optional
Retain original materials within primitive and fill outside.
Defaults to False.
periodic : Boolean, optional
Repeat primitive over boundaries. Defaults to True.
"""
# normalized 'radius' and center
r = np.array(dimension)/self.grid/2.0 if np.array(dimension).dtype in np.sctypes['int'] else \
np.array(dimension)/self.size/2.0
c = (np.array(center) + .5)/self.grid if np.array(center).dtype in np.sctypes['int'] else \
(np.array(center) - self.origin)/self.size
coords = grid_filters.cell_coord0(self.grid,np.ones(3)) \
- ((np.ones(3)-(1./self.grid if np.array(center).dtype in np.sctypes['int'] else 0))*0.5 if periodic else c) # periodic center is always at CoG
coords_rot = R.broadcast_to(tuple(self.grid))@coords
with np.errstate(all='ignore'):
mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) > 1.0,True,False)
if periodic: # translate back to center
mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))
fill_ = np.full_like(self.material,np.nanmax(self.material)+1 if fill is None else fill)
return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','add_primitive')],
)
def mirror(self,directions,reflect=False):
"""
Mirror geometry along given directions.
Parameters
----------
directions : iterable containing str
Direction(s) along which the geometry is mirrored.
Valid entries are 'x', 'y', 'z'.
reflect : bool, optional
Reflect (include) outermost layers. Defaults to False.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
limits = [None,None] if reflect else [-2,0]
mat = self.material.copy()
if 'x' in directions:
mat = np.concatenate([mat,mat[limits[0]:limits[1]:-1,:,:]],0)
if 'y' in directions:
mat = np.concatenate([mat,mat[:,limits[0]:limits[1]:-1,:]],1)
if 'z' in directions:
mat = np.concatenate([mat,mat[:,:,limits[0]:limits[1]:-1]],2)
return Geom(material = mat,
size = self.size/self.grid*np.asarray(mat.shape),
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','mirror')],
)
def flip(self,directions):
"""
Flip geometry along given directions.
Parameters
----------
directions : iterable containing str
Direction(s) along which the geometry is flipped.
Valid entries are 'x', 'y', 'z'.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
return Geom(material = mat,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','flip')],
)
def scale(self,grid,periodic=True):
"""
Scale geometry to new grid.
Parameters
----------
grid : numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
periodic : Boolean, optional
Assume geometry to be periodic. Defaults to True.
"""
return Geom(material = ndimage.interpolation.zoom(
self.material,
grid/self.grid,
output=self.material.dtype,
order=0,
mode=('wrap' if periodic else 'nearest'),
prefilter=False
),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','scale')],
)
def clean(self,stencil=3,selection=None,periodic=True):
"""
Smooth geometry by selecting most frequent material index within given stencil at each location.
Parameters
----------
stencil : int, optional
Size of smoothing stencil.
selection : list, optional
Field values that can be altered. Defaults to all.
periodic : Boolean, optional
Assume geometry to be periodic. Defaults to True.
"""
def mostFrequent(arr,selection=None):
me = arr[arr.size//2]
if selection is None or me in selection:
unique, inverse = np.unique(arr, return_inverse=True)
return unique[np.argmax(np.bincount(inverse))]
else:
return me
return Geom(material = ndimage.filters.generic_filter(
self.material,
mostFrequent,
size=(stencil if selection is None else stencil//2*2+1,)*3,
mode=('wrap' if periodic else 'nearest'),
extra_keywords=dict(selection=selection),
).astype(self.material.dtype),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','clean')],
)
def renumber(self):
"""Renumber sorted material indices to 1,...,N."""
renumbered = np.empty(self.grid,dtype=self.material.dtype)
for i, oldID in enumerate(np.unique(self.material)):
renumbered = np.where(self.material == oldID, i+1, renumbered)
return Geom(material = renumbered,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','renumber')],
)
def rotate(self,R,fill=None):
"""
Rotate geometry (pad if required).
Parameters
----------
R : damask.Rotation
Rotation to apply to the geometry.
fill : int or float, optional
Material index to fill the corners. Defaults to material.max() + 1.
"""
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
Eulers = R.as_Eulers(degrees=True)
material_in = self.material.copy()
# 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
for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]):
material_out = ndimage.rotate(material_in,angle,axes,order=0,
prefilter=False,output=dtype,cval=fill)
if np.prod(material_in.shape) == np.prod(material_out.shape):
# avoid scipy interpolation errors for rotations close to multiples of 90°
material_in = np.rot90(material_in,k=np.rint(angle/90.).astype(int),axes=axes)
else:
material_in = material_out
origin = self.origin-(np.asarray(material_in.shape)-self.grid)*.5 * self.size/self.grid
return Geom(material = material_in,
size = self.size/self.grid*np.asarray(material_in.shape),
origin = origin,
comments = self.comments+[util.execution_stamp('Geom','rotate')],
)
def canvas(self,grid=None,offset=None,fill=None):
"""
Crop or enlarge/pad geometry.
Parameters
----------
grid : numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
offset : numpy.ndarray of shape (3)
Offset (measured in grid points) from old to new geometry [0,0,0].
fill : int or float, optional
Material index to fill the background. Defaults to material.max() + 1.
"""
if offset is None: offset = 0
if fill is None: fill = np.nanmax(self.material) + 1
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)
LL = np.clip( offset, 0,np.minimum(self.grid, grid+offset))
UR = np.clip( offset+grid, 0,np.minimum(self.grid, 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))
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(material = canvas,
size = self.size/self.grid*np.asarray(canvas.shape),
origin = self.origin+offset*self.size/self.grid,
comments = self.comments+[util.execution_stamp('Geom','canvas')],
)
def substitute(self,from_material,to_material):
"""
Substitute material indices.
Parameters
----------
from_material : iterable of ints
Material indices to be substituted.
to_material : iterable of ints
New material indices.
"""
substituted = self.material.copy()
for from_ms,to_ms in zip(from_material,to_material):
substituted[self.material==from_ms] = to_ms
return Geom(material = substituted,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','substitute')],
)
def vicinity_offset(self,vicinity=1,offset=None,trigger=[],periodic=True):
"""
Offset material index of points in the vicinity of xxx.
Different from themselves (or listed as triggers) within a given (cubic) vicinity,
i.e. within the region close to a grain/phase boundary.
ToDo: use include/exclude as in seeds.from_geom
Parameters
----------
vicinity : int, optional
Voxel distance checked for presence of other materials.
Defaults to 1.
offset : int, optional
Offset (positive or negative) to tag material indices,
defaults to material.max() + 1.
trigger : list of ints, optional
List of material indices that trigger a change.
Defaults to [], meaning that any different neighbor triggers a change.
periodic : Boolean, optional
Assume geometry to be periodic. Defaults to True.
"""
def tainted_neighborhood(stencil,trigger):
me = stencil[stencil.shape[0]//2]
if len(trigger) == 0:
return np.any(stencil != me)
if me in trigger:
trigger = set(trigger)
trigger.remove(me)
trigger = list(trigger)
return np.any(np.in1d(stencil,np.array(trigger)))
offset_ = np.nanmax(self.material) if offset is None else offset
mask = ndimage.filters.generic_filter(self.material,
tainted_neighborhood,
size=1+2*vicinity,
mode='wrap' if periodic else 'nearest',
extra_keywords={'trigger':trigger})
return Geom(material = np.where(mask, self.material + offset_,self.material),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')],
)