now using ASCIItable object (much faster file writing).
(temporarily) switched back to rev2496 logic... changed to diffusion distance as command line argument (instead of taking sqrt(time))
This commit is contained in:
Philip Eisenlohr
parent
dd3d53e238
commit
1f891c544d
|
|
@ -2,6 +2,7 @@
|
|||
# -*- coding: UTF-8 no BOM -*-
|
||||
|
||||
import os,sys,string,re,math,numpy
|
||||
import damask
|
||||
from optparse import OptionParser, OptionGroup, Option, SUPPRESS_HELP
|
||||
from scipy import ndimage
|
||||
|
||||
|
|
@ -28,37 +29,37 @@ class extendedOption(Option):
|
|||
# MAIN
|
||||
#--------------------------------------------------------------------------------------------------
|
||||
|
||||
synonyms = {
|
||||
'grid': ['resolution'],
|
||||
'size': ['dimension'],
|
||||
}
|
||||
identifiers = {
|
||||
'grid': ['a','b','c'],
|
||||
'size': ['x','y','z'],
|
||||
'origin': ['x','y','z'],
|
||||
}
|
||||
mappings = {
|
||||
'grid': lambda x: int(x),
|
||||
'size': lambda x: float(x),
|
||||
'origin': lambda x: float(x),
|
||||
'homogenization': lambda x: int(x),
|
||||
'grid': lambda x: int(x),
|
||||
'size': lambda x: float(x),
|
||||
'origin': lambda x: float(x),
|
||||
'homogenization': lambda x: int(x),
|
||||
'microstructures': lambda x: int(x),
|
||||
}
|
||||
|
||||
parser = OptionParser(option_class=extendedOption, usage='%prog options [file[s]]', description = """
|
||||
Smoothens out interface roughness by simulated curvature flow.
|
||||
This is achieved by the diffusion of each initially sharply bounded grain volume within the periodic domain for a given time,
|
||||
i.e. up to a diffusion distance of sqrt(t) voxels.
|
||||
The final geometry is assembled by selecting at each voxel that grain index for which the concentration is largest.
|
||||
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.
|
||||
""" + string.replace('$Id$','\n','\\n')
|
||||
)
|
||||
|
||||
parser.add_option('-t', '--time', dest='t', type='int', \
|
||||
help='time for curvature flow [%default]')
|
||||
parser.add_option('-N', '--smooth', dest='N', type='int', \
|
||||
help='number of steps for curvature flow [%default]')
|
||||
parser.add_option('-d', '--distance', dest='d', type='float', \
|
||||
help='diffusion distance in voxels [%default]', metavar='float')
|
||||
parser.add_option('-b', '--black', dest='black', action='extend', type='string', \
|
||||
help='indices of stationary microstructures', metavar='<LIST>')
|
||||
parser.add_option('-2', '--twodimensional', dest='twoD', action='store_true', \
|
||||
help='output geom file with two-dimensional data arrangement [%default]')
|
||||
|
||||
parser.set_defaults(t = 1)
|
||||
parser.set_defaults(N = 1)
|
||||
parser.set_defaults(black = [])
|
||||
parser.set_defaults(twoD = False)
|
||||
|
||||
|
|
@ -83,36 +84,31 @@ else:
|
|||
'croak':sys.stdout,
|
||||
})
|
||||
|
||||
#--- loop over input files ------------------------------------------------------------------------
|
||||
#--- loop over input files ------------------------------------------------------------------------
|
||||
for file in files:
|
||||
if file['name'] != 'STDIN': file['croak'].write(file['name']+'\n')
|
||||
|
||||
firstline = file['input'].readline()
|
||||
m = re.search('(\d+)\s*head', firstline.lower())
|
||||
if m:
|
||||
headerlines = int(m.group(1))
|
||||
headers = [file['input'].readline() for i in range(headerlines)]
|
||||
else:
|
||||
headerlines = 1
|
||||
headers = firstline
|
||||
theTable = damask.ASCIItable(file['input'],file['output'],labels=False)
|
||||
theTable.head_read()
|
||||
|
||||
content = file['input'].readlines()
|
||||
file['input'].close()
|
||||
|
||||
#--- interprete header ----------------------------------------------------------------------------
|
||||
#--- interpret header ----------------------------------------------------------------------------
|
||||
info = {
|
||||
'grid': numpy.zeros(3,'i'),
|
||||
'size': numpy.zeros(3,'d'),
|
||||
'origin': numpy.zeros(3,'d'),
|
||||
'microstructures': 0,
|
||||
'grid': numpy.zeros(3,'i'),
|
||||
'size': numpy.zeros(3,'d'),
|
||||
'origin': numpy.zeros(3,'d'),
|
||||
'homogenization': 0,
|
||||
'microstructures': 0,
|
||||
}
|
||||
newInfo = {
|
||||
'microstructures': 0,
|
||||
}
|
||||
extra_header = []
|
||||
|
||||
new_header = []
|
||||
for header in headers:
|
||||
for header in theTable.info:
|
||||
headitems = map(str.lower,header.split())
|
||||
if headitems[0] == 'resolution': headitems[0] = 'grid'
|
||||
if headitems[0] == 'dimension': headitems[0] = 'size'
|
||||
if len(headitems) == 0: continue
|
||||
for synonym,alternatives in synonyms.iteritems():
|
||||
if headitems[0] in alternatives: headitems[0] = synonym
|
||||
if headitems[0] in mappings.keys():
|
||||
if headitems[0] in identifiers.keys():
|
||||
for i in xrange(len(identifiers[headitems[0]])):
|
||||
|
|
@ -120,6 +116,8 @@ for file in files:
|
|||
mappings[headitems[0]](headitems[headitems.index(identifiers[headitems[0]][i])+1])
|
||||
else:
|
||||
info[headitems[0]] = mappings[headitems[0]](headitems[1])
|
||||
else:
|
||||
extra_header.append(header)
|
||||
|
||||
file['croak'].write('grid a b c: %s\n'%(' x '.join(map(str,info['grid']))) + \
|
||||
'size x y z: %s\n'%(' x '.join(map(str,info['size']))) + \
|
||||
|
|
@ -129,100 +127,69 @@ for file in files:
|
|||
|
||||
if numpy.any(info['grid'] < 1):
|
||||
file['croak'].write('invalid grid a b c.\n')
|
||||
sys.exit()
|
||||
continue
|
||||
if numpy.any(info['size'] <= 0.0):
|
||||
file['croak'].write('invalid size x y z.\n')
|
||||
sys.exit()
|
||||
continue
|
||||
|
||||
#--- read data ------------------------------------------------------------------------------------
|
||||
microstructure = numpy.zeros(info['grid'],'i')
|
||||
#--- read data ------------------------------------------------------------------------------------
|
||||
microstructure = numpy.zeros(info['grid'].prod(),'i')
|
||||
i = 0
|
||||
for line in content:
|
||||
items = line.split()
|
||||
theTable.data_rewind()
|
||||
while theTable.data_read():
|
||||
items = theTable.data
|
||||
if len(items) > 2:
|
||||
if items[1].lower() == 'of': items = [int(items[2])]*int(items[0])
|
||||
elif items[1].lower() == 'to': items = xrange(int(items[0]),1+int(items[2]))
|
||||
else: items = map(int,items)
|
||||
else: items = map(int,items)
|
||||
else: items = map(int,items)
|
||||
else: items = map(int,items)
|
||||
|
||||
for item in items:
|
||||
microstructure[i%info['grid'][0],
|
||||
(i/info['grid'][0])%info['grid'][1],
|
||||
i/info['grid'][0] /info['grid'][1]] = item
|
||||
i += 1
|
||||
s = len(items)
|
||||
microstructure[i:i+s] = items
|
||||
i += s
|
||||
|
||||
#--- do work ------------------------------------------------------------------------------------
|
||||
microstructure = microstructure.reshape(info['grid'],order='F')
|
||||
|
||||
#--- initialize helper data -----------------------------------------------------------------------
|
||||
|
||||
diffusionWindow = int(math.ceil(4.0*numpy.sqrt(options.t)))
|
||||
for smoothIter in xrange(options.N):
|
||||
extendedMicro = numpy.zeros(2*diffusionWindow+info['grid']).astype(int)
|
||||
extendedMicro[:info['grid'][0],:info['grid'][1],:info['grid'][2]] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=diffusionWindow,axis=0),\
|
||||
shift=diffusionWindow,axis=1),shift=diffusionWindow,axis=2)
|
||||
extendedMicro[-info['grid'][0]:,:info['grid'][1],:info['grid'][2]] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=-diffusionWindow,axis=0),\
|
||||
shift=diffusionWindow,axis=1),shift=diffusionWindow,axis=2)
|
||||
extendedMicro[:info['grid'][0],-info['grid'][1]:,:info['grid'][2]] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=diffusionWindow,axis=0),\
|
||||
shift=-diffusionWindow,axis=1),shift=diffusionWindow,axis=2)
|
||||
extendedMicro[-info['grid'][0]:,-info['grid'][1]:,:info['grid'][2]] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=-diffusionWindow,axis=0),\
|
||||
shift=-diffusionWindow,axis=1),shift=diffusionWindow,axis=2)
|
||||
extendedMicro[:info['grid'][0],:info['grid'][1],-info['grid'][2]:] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=diffusionWindow,axis=0),\
|
||||
shift=diffusionWindow,axis=1),shift=-diffusionWindow,axis=2)
|
||||
extendedMicro[-info['grid'][0]:,:info['grid'][1],-info['grid'][2]:] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=-diffusionWindow,axis=0),\
|
||||
shift=diffusionWindow,axis=1),shift=-diffusionWindow,axis=2)
|
||||
extendedMicro[:info['grid'][0],-info['grid'][1]:,-info['grid'][2]:] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=diffusionWindow,axis=0),\
|
||||
shift=-diffusionWindow,axis=1),shift=-diffusionWindow,axis=2)
|
||||
extendedMicro[-info['grid'][0]:,-info['grid'][1]:,-info['grid'][2]:] = \
|
||||
numpy.roll(numpy.roll(numpy.roll(microstructure,shift=-diffusionWindow,axis=0),\
|
||||
shift=-diffusionWindow,axis=1),shift=-diffusionWindow,axis=2)
|
||||
winner = numpy.zeros(extendedMicro.shape).astype(int)
|
||||
winner[diffusionWindow:-diffusionWindow,diffusionWindow:-diffusionWindow,diffusionWindow:-diffusionWindow] =\
|
||||
numpy.where(numpy.reshape(numpy.in1d(microstructure,options.black),microstructure.shape),\
|
||||
microstructure,0)
|
||||
diffusedMax = numpy.zeros(extendedMicro.shape)
|
||||
boundingSlice = ndimage.measurements.find_objects(microstructure)
|
||||
microList = set(numpy.unique(microstructure)).difference(set(options.black).union(set([0])))
|
||||
winner = numpy.zeros(info['grid'],'i')
|
||||
diffusedMax = numpy.zeros(info['grid'])
|
||||
|
||||
#--- diffuse each grain separately ----------------------------------------------------------------
|
||||
for theGrain in xrange(1,1+numpy.amax(microstructure)):
|
||||
diffused = ndimage.filters.gaussian_filter((microstructure == theGrain).astype(float),\
|
||||
{True:0.0,False:options.d}[theGrain in options.black],\
|
||||
mode='wrap')
|
||||
winner = numpy.where(diffused > diffusedMax, theGrain, winner)
|
||||
diffusedMax = numpy.where(diffused > diffusedMax, diffused, diffusedMax)
|
||||
|
||||
for grain in microList:
|
||||
xMin = boundingSlice[grain-1][0].start; xMax = boundingSlice[grain-1][0].stop + 2*diffusionWindow
|
||||
yMin = boundingSlice[grain-1][1].start; yMax = boundingSlice[grain-1][1].stop + 2*diffusionWindow
|
||||
zMin = boundingSlice[grain-1][2].start; zMax = boundingSlice[grain-1][2].stop + 2*diffusionWindow
|
||||
diffused = ndimage.filters.gaussian_filter((extendedMicro[xMin:xMax,yMin:yMax,zMin:zMax] == grain).astype(float),\
|
||||
numpy.sqrt(options.t))
|
||||
isMax = diffused > diffusedMax[xMin:xMax,yMin:yMax,zMin:zMax]
|
||||
winner[xMin:xMax,yMin:yMax,zMin:zMax][isMax] = grain
|
||||
diffusedMax[xMin:xMax,yMin:yMax,zMin:zMax] = numpy.where(isMax,diffused,diffusedMax[xMin:xMax,yMin:yMax,zMin:zMax])
|
||||
|
||||
microstructure = winner[diffusionWindow:-diffusionWindow,diffusionWindow:-diffusionWindow,diffusionWindow:-diffusionWindow]
|
||||
|
||||
# --- assemble header -----------------------------------------------------------------------------
|
||||
formatwidth = int(math.floor(math.log10(microstructure.max())+1))
|
||||
newInfo['microstructures'] = winner.max()
|
||||
|
||||
new_header.append('$Id$\n')
|
||||
new_header.append("grid\ta %i\tb %i\tc %i\n"%(
|
||||
info['grid'][0],info['grid'][1],info['grid'][2]))
|
||||
new_header.append("size\tx %f\ty %f\tz %f\n"%(
|
||||
info['size'][0],info['size'][1],info['size'][2]))
|
||||
new_header.append("origin\tx %f\ty %f\tz %f\n"%(
|
||||
info['origin'][0],info['origin'][1],info['origin'][2]))
|
||||
new_header.append("homogenization\t%i\n"%info['homogenization'])
|
||||
new_header.append("microstructures\t%i\n"%info['microstructures'])
|
||||
file['output'].write('%i\theader\n'%(len(new_header))+''.join(new_header))
|
||||
#--- report ---------------------------------------------------------------------------------------
|
||||
if (newInfo['microstructures'] != info['microstructures']):
|
||||
file['croak'].write('--> microstructures: %i\n'%newInfo['microstructures'])
|
||||
|
||||
#--- write header ---------------------------------------------------------------------------------
|
||||
theTable.labels_clear()
|
||||
theTable.info_clear()
|
||||
theTable.info_append(extra_header+[
|
||||
"$Id$",
|
||||
"grid\ta %i\tb %i\tc %i"%(info['grid'][0],info['grid'][1],info['grid'][2],),
|
||||
"size\tx %f\ty %f\tz %f"%(info['size'][0],info['size'][1],info['size'][2],),
|
||||
"origin\tx %f\ty %f\tz %f"%(info['origin'][0],info['origin'][1],info['origin'][2],),
|
||||
"homogenization\t%i"%info['homogenization'],
|
||||
"microstructures\t%i"%(newInfo['microstructures']),
|
||||
])
|
||||
theTable.head_write()
|
||||
theTable.output_flush()
|
||||
|
||||
# --- write microstructure information ------------------------------------------------------------
|
||||
for z in xrange(info['grid'][2]):
|
||||
for y in xrange(info['grid'][1]):
|
||||
file['output'].write({True:' ',False:'\n'}[options.twoD].join(map(lambda x: \
|
||||
('%%%ii'%formatwidth)%x, microstructure[:,y,z])) + '\n')
|
||||
|
||||
formatwidth = int(math.floor(math.log10(winner.max())+1))
|
||||
theTable.data = winner.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose()
|
||||
theTable.data_writeArray('%%%ii'%(formatwidth))
|
||||
|
||||
#--- output finalization --------------------------------------------------------------------------
|
||||
if file['name'] != 'STDIN':
|
||||
file['input'].close()
|
||||
file['output'].close()
|
||||
os.rename(file['name']+'_tmp',file['name'])
|
||||
|
|
|
|||
Loading…
Reference in New Issue