From 1f891c544d467a750b70f94ff756285f86b660fb Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 30 Jun 2013 13:51:21 +0000 Subject: [PATCH] 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)) --- processing/pre/geom_stretchInterfaces.py | 191 ++++++++++------------- 1 file changed, 79 insertions(+), 112 deletions(-) diff --git a/processing/pre/geom_stretchInterfaces.py b/processing/pre/geom_stretchInterfaces.py index 348eb97ba..3d4f2ea10 100755 --- a/processing/pre/geom_stretchInterfaces.py +++ b/processing/pre/geom_stretchInterfaces.py @@ -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='') -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'])