DAMASK_EICMD/processing/pre/geom_fromVoronoiTessellatio...

342 lines
14 KiB
Python
Raw Normal View History

#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
2016-03-02 15:59:07 +05:30
import os,sys,math
import numpy as np
import multiprocessing
from optparse import OptionParser,OptionGroup
from scipy import spatial
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
2014-10-09 22:33:06 +05:30
def meshgrid2(*arrs):
2016-03-02 15:59:07 +05:30
"""code inspired by http://stackoverflow.com/questions/1827489/numpy-meshgrid-in-3d"""
arrs = tuple(reversed(arrs))
arrs = tuple(arrs)
lens = np.array(map(len, arrs))
dim = len(arrs)
ans = []
for i, arr in enumerate(arrs):
slc = np.ones(dim,'i')
slc[i] = lens[i]
arr2 = np.asarray(arr).reshape(slc)
for j, sz in enumerate(lens):
if j != i:
arr2 = arr2.repeat(sz, axis=j)
ans.insert(0,arr2)
return tuple(ans)
def findClosestSeed(fargs):
2016-07-15 19:18:53 +05:30
point, seeds, myWeights = fargs
tmp = np.repeat(point.reshape(3,1), len(seeds), axis=1).T
2016-07-15 19:18:53 +05:30
dist = np.sum((tmp - seeds)**2,axis=1) -myWeights
return np.argmin(dist) # seed point closest to point
def laguerreTessellation(undeformed, coords, weights, grains, nonperiodic = False, cpus = 2):
copies = \
np.array([
[ 0, 0, 0 ],
]).astype(float) if nonperiodic else \
2016-03-02 15:59:07 +05:30
np.array([
[ -1,-1,-1 ],
[ 0,-1,-1 ],
[ 1,-1,-1 ],
[ -1, 0,-1 ],
[ 0, 0,-1 ],
[ 1, 0,-1 ],
[ -1, 1,-1 ],
[ 0, 1,-1 ],
[ 1, 1,-1 ],
[ -1,-1, 0 ],
[ 0,-1, 0 ],
[ 1,-1, 0 ],
[ -1, 0, 0 ],
[ 0, 0, 0 ],
[ 1, 0, 0 ],
[ -1, 1, 0 ],
[ 0, 1, 0 ],
[ 1, 1, 0 ],
[ -1,-1, 1 ],
[ 0,-1, 1 ],
[ 1,-1, 1 ],
[ -1, 0, 1 ],
[ 0, 0, 1 ],
[ 1, 0, 1 ],
[ -1, 1, 1 ],
[ 0, 1, 1 ],
[ 1, 1, 1 ],
]).astype(float)*info['size']
repeatweights = np.tile(weights,len(copies)).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
for i,vec in enumerate(copies): # periodic copies of seed points ...
try: seeds = np.append(seeds, coords+vec, axis=0) # ... (1+a,2+a,3+a,...,1+z,2+z,3+z)
except NameError: seeds = coords+vec
if (repeatweights == 0.0).all(): # standard Voronoi (no weights, KD tree)
myKDTree = spatial.cKDTree(seeds)
devNull,closestSeeds = myKDTree.query(undeformed)
else:
damask.util.croak('...using {} cpu{}'.format(options.cpus, 's' if options.cpus > 1 else ''))
2016-07-15 18:22:48 +05:30
arguments = [[arg,seeds,repeatweights] for arg in list(undeformed)]
if cpus > 1: # use multithreading
pool = multiprocessing.Pool(processes = cpus) # initialize workers
result = pool.map_async(findClosestSeed, arguments) # evaluate function in parallel
pool.close()
pool.join()
closestSeeds = np.array(result.get()).flatten()
else:
closestSeeds = np.zeros(len(arguments),dtype='i')
for i,arg in enumerate(arguments):
closestSeeds[i] = findClosestSeed(arg)
2016-03-02 15:59:07 +05:30
# closestSeed is modulo number of original seed points (i.e. excluding periodic copies)
return grains[closestSeeds%coords.shape[0]]
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
2014-10-09 22:33:06 +05:30
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [seedfile(s)]', description = """
Generate geometry description and material configuration by tessellation of given seeds file.
""", version = scriptID)
group = OptionGroup(parser, "Tessellation","")
group.add_option('-l',
'--laguerre',
dest = 'laguerre',
action = 'store_true',
help = 'use Laguerre (weighted Voronoi) tessellation')
group.add_option('--cpus',
dest = 'cpus',
type = 'int', metavar = 'int',
help = 'number of parallel processes to use for Laguerre tessellation [%default]')
group.add_option('--nonperiodic',
dest = 'nonperiodic',
action = 'store_true',
help = 'nonperiodic tessellation')
parser.add_option_group(group)
group = OptionGroup(parser, "Geometry","")
group.add_option('-g',
'--grid',
dest = 'grid',
type = 'int', nargs = 3, metavar = ' '.join(['int']*3),
help = 'a,b,c grid of hexahedral box')
group.add_option('-s',
'--size',
dest = 'size',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'x,y,z size of hexahedral box')
group.add_option('-o',
'--origin',
dest = 'origin',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'origin of grid')
parser.add_option_group(group)
group = OptionGroup(parser, "Seeds","")
group.add_option('-p',
'--pos', '--seedposition',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
group.add_option('-w',
'--weight',
dest = 'weight',
type = 'string', metavar = 'string',
help = 'label of weights [%default]')
group.add_option('-m',
'--microstructure',
dest = 'microstructure',
type = 'string', metavar = 'string',
help = 'label of microstructures [%default]')
group.add_option('-e',
'--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'label of Euler angles [%default]')
group.add_option('--axes',
dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame')
parser.add_option_group(group)
group = OptionGroup(parser, "Configuration","")
group.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
group.add_option('--crystallite',
dest = 'crystallite',
type = 'int', metavar = 'int',
help = 'crystallite index to be used [%default]')
group.add_option('--phase',
dest = 'phase',
type = 'int', metavar = 'int',
help = 'phase index to be used [%default]')
parser.add_option_group(group)
parser.set_defaults(pos = 'pos',
weight = 'weight',
microstructure = 'microstructure',
eulers = 'euler',
homogenization = 1,
crystallite = 1,
phase = 1,
cpus = 2,
laguerre = False,
nonperiodic = False,
)
(options,filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
outname = os.path.splitext(name)[0]+'.geom' if name else name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# --- read header ----------------------------------------------------------------------------
table.head_read()
info,extra_header = table.head_getGeom()
2016-03-02 15:59:07 +05:30
if options.grid is not None: info['grid'] = options.grid
if options.size is not None: info['size'] = options.size
if options.origin is not None: info['origin'] = options.origin
# ------------------------------------------ sanity checks ---------------------------------------
remarks = []
errors = []
labels = []
hasGrains = table.label_dimension(options.microstructure) == 1
hasEulers = table.label_dimension(options.eulers) == 3
hasWeights = table.label_dimension(options.weight) == 1 and options.laguerre
if np.any(info['grid'] < 1): errors.append('invalid grid a b c.')
if np.any(info['size'] <= 0.0) \
and np.all(info['grid'] < 1): errors.append('invalid size x y z.')
else:
for i in xrange(3):
if info['size'][i] <= 0.0: # any invalid size?
info['size'][i] = float(info['grid'][i])/max(info['grid']) # normalize to grid
remarks.append('rescaling size {} to {}...'.format({0:'x',1:'y',2:'z'}[i],info['size'][i]))
if table.label_dimension(options.pos) != 3:
errors.append('seed positions "{}" have dimension {}.'.format(options.pos,
table.label_dimension(options.pos)))
else:
labels += [options.pos]
if not hasEulers: remarks.append('missing seed orientations...')
else: labels += [options.eulers]
if not hasGrains: remarks.append('missing seed microstructure indices...')
else: labels += [options.microstructure]
if options.laguerre and not hasWeights: remarks.append('missing seed weights...')
else: labels += [options.weight]
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss=True)
continue
# ------------------------------------------ read seeds ---------------------------------------
table.data_readArray(labels)
coords = table.data[:,table.label_indexrange(options.pos)] * info['size']
eulers = table.data[:,table.label_indexrange(options.eulers)] if hasEulers \
else np.zeros(3*len(coords))
grains = table.data[:,table.label_indexrange(options.microstructure)].astype('i') if hasGrains \
else 1+np.arange(len(coords))
weights = table.data[:,table.label_indexrange(options.weight)] if hasWeights \
else np.zeros(len(coords))
grainIDs = np.unique(grains).astype('i')
NgrainIDs = len(grainIDs)
# --- tessellate microstructure ------------------------------------------------------------
x = (np.arange(info['grid'][0])+0.5)*info['size'][0]/info['grid'][0]
y = (np.arange(info['grid'][1])+0.5)*info['size'][1]/info['grid'][1]
z = (np.arange(info['grid'][2])+0.5)*info['size'][2]/info['grid'][2]
damask.util.croak('tessellating...')
grid = np.vstack(meshgrid2(x, y, z)).reshape(3,-1).T
indices = laguerreTessellation(grid, coords, weights, grains, options.nonperiodic, options.cpus)
# --- write header ------------------------------------------------------------------------
usedGrainIDs = np.intersect1d(grainIDs,indices)
info['microstructures'] = len(usedGrainIDs)
if info['homogenization'] == 0: info['homogenization'] = options.homogenization
damask.util.report_geom(info,['grid','size','origin','homogenization',])
damask.util.croak(['microstructures: {}{}'.format(info['microstructures'],
(' out of {}'.format(NgrainIDs) if NgrainIDs != info['microstructures'] else '')),
])
config_header = []
formatwidth = 1+int(math.log10(NgrainIDs))
config_header += ['<microstructure>']
for i,ID in enumerate(grainIDs):
config_header += ['[Grain{}]'.format(str(ID).zfill(formatwidth)),
'crystallite {}'.format(options.crystallite),
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(options.phase,str(ID).rjust(formatwidth)),
]
if hasEulers:
config_header += ['<texture>']
for ID in grainIDs:
eulerID = np.nonzero(grains == ID)[0][0] # find first occurrence of this grain id
config_header += ['[Grain{}]'.format(str(ID).zfill(formatwidth)),
'(gauss)\tphi1 {:g}\tPhi {:g}\tphi2 {:g}\tscatter 0.0\tfraction 1.0'.format(*eulers[eulerID])
]
if options.axes is not None: config_header.append('axes\t{} {} {}'.format(*options.axes))
table.labels_clear()
table.info_clear()
table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {}\tb {}\tc {}".format(*info['grid']),
"size\tx {}\ty {}\tz {}".format(*info['size']),
"origin\tx {}\ty {}\tz {}".format(*info['origin']),
"homogenization\t{}".format(info['homogenization']),
"microstructures\t{}".format(info['microstructures']),
config_header,
])
table.head_write()
# --- write microstructure information ------------------------------------------------------------
table.data = indices.reshape(info['grid'][1]*info['grid'][2],info['grid'][0])
table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
#--- output finalization --------------------------------------------------------------------------
table.close()