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

fixed error in file handling

This commit is contained in:
Martin Diehl 2014-06-18 09:00:57 +00:00
parent 69857cc608
commit 1ff38b98af
1 changed files with 35 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
processing/pre/geom_grainGrowth.py
View File

@ -1,34 +1,16 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys,string,re,math,numpy,itertools import os,sys,string,re,math,itertools
import damask import numpy as np
from optparse import OptionParser, OptionGroup, Option, SUPPRESS_HELP from optparse import OptionParser
from scipy import ndimage from scipy import ndimage
from multiprocessing import Pool from multiprocessing import Pool
import damask
scriptID = '$Id$' scriptID = '$Id$'
scriptName = scriptID.split()[1] scriptName = scriptID.split()[1]
#--------------------------------------------------------------------------------------------------
class extendedOption(Option):
#--------------------------------------------------------------------------------------------------
# used for definition of new option parser action 'extend', which enables to take multiple option arguments
# taken from online tutorial http://docs.python.org/library/optparse.html
ACTIONS = Option.ACTIONS + ("extend",)
STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# MAIN # MAIN
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
@ -49,7 +31,7 @@ mappings = {
'microstructures': lambda x: int(x), 'microstructures': lambda x: int(x),
} }
parser = OptionParser(option_class=extendedOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Smoothens out interface roughness by simulated curvature flow. 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 This is achieved by the diffusion of each initially sharply bounded grain volume within the periodic domain
up to a given distance 'd' voxels. up to a given distance 'd' voxels.
@ -102,9 +84,9 @@ for file in files:
#--- interpret header ---------------------------------------------------------------------------- #--- interpret header ----------------------------------------------------------------------------
info = { info = {
'grid': numpy.zeros(3,'i'), 'grid': np.zeros(3,'i'),
'size': numpy.zeros(3,'d'), 'size': np.zeros(3,'d'),
'origin': numpy.zeros(3,'d'), 'origin': np.zeros(3,'d'),
'homogenization': 0, 'homogenization': 0,
'microstructures': 0, 'microstructures': 0,
} }
@ -134,15 +116,15 @@ for file in files:
'homogenization: %i\n'%info['homogenization'] + \ 'homogenization: %i\n'%info['homogenization'] + \
'microstructures: %i\n'%info['microstructures']) 'microstructures: %i\n'%info['microstructures'])
if numpy.any(info['grid'] < 1): if np.any(info['grid'] < 1):
file['croak'].write('invalid grid a b c.\n') file['croak'].write('invalid grid a b c.\n')
continue continue
if numpy.any(info['size'] <= 0.0): if np.any(info['size'] <= 0.0):
file['croak'].write('invalid size x y z.\n') file['croak'].write('invalid size x y z.\n')
continue continue
#--- read data ------------------------------------------------------------------------------------ #--- read data ------------------------------------------------------------------------------------
microstructure = numpy.zeros(numpy.prod([2 if i == 1 else i for i in info['grid']]),'i') # 2D structures do not work microstructure = np.zeros(np.prod([2 if i == 1 else i for i in info['grid']]),'i') # 2D structures do not work
i = 0 i = 0
while theTable.data_read(): # read next data line of ASCII table while theTable.data_read(): # read next data line of ASCII table
@ -159,47 +141,47 @@ for file in files:
#--- reshape, if 2D make copy --------------------------------------------------------------------- #--- reshape, if 2D make copy ---------------------------------------------------------------------
microstructure = numpy.tile(microstructure.reshape(info['grid'],order='F'), microstructure = np.tile(microstructure.reshape(info['grid'],order='F'),
numpy.where(info['grid'] == 1, 2,1)) # make one copy along dimensions with grid == 1 np.where(info['grid'] == 1, 2,1)) # make one copy along dimensions with grid == 1
grid = numpy.array(microstructure.shape) grid = np.array(microstructure.shape)
#--- initialize support data ----------------------------------------------------------------------- #--- initialize support data -----------------------------------------------------------------------
periodic_microstructure = numpy.tile(microstructure,(3,3,3))[grid[0]/2:-grid[0]/2, periodic_microstructure = np.tile(microstructure,(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 microstructure grid[2]/2:-grid[2]/2] # periodically extend the microstructure
microstructure_original = numpy.copy(microstructure) # store a copy the initial microstructure to find locations of immutable indices microstructure_original = np.copy(microstructure) # store a copy the initial microstructure to find locations of immutable indices
X,Y,Z = numpy.mgrid[0:grid[0],0:grid[1],0:grid[2]] X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]]
gauss = numpy.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d))/math.pow(2.0*numpy.pi*options.d*options.d,1.5) gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d))/math.pow(2.0*np.pi*options.d*options.d,1.5)
gauss[:,:,grid[2]/2::] = gauss[:,:,round(grid[2]/2.)-1::-1] # trying to cope with uneven (odd) grid size gauss[:,:,grid[2]/2::] = gauss[:,:,round(grid[2]/2.)-1::-1] # trying to cope with uneven (odd) grid size
gauss[:,grid[1]/2::,:] = gauss[:,round(grid[1]/2.)-1::-1,:] gauss[:,grid[1]/2::,:] = gauss[:,round(grid[1]/2.)-1::-1,:]
gauss[grid[0]/2::,:,:] = gauss[round(grid[0]/2.)-1::-1,:,:] gauss[grid[0]/2::,:,:] = gauss[round(grid[0]/2.)-1::-1,:,:]
gauss = numpy.fft.rfftn(gauss) gauss = np.fft.rfftn(gauss)
interfacialEnergy = lambda A,B: (A*B != 0)*(A != B)*1.0 interfacialEnergy = lambda A,B: (A*B != 0)*(A != B)*1.0
struc = ndimage.generate_binary_structure(3,1) # 3D von Neumann neighborhood struc = ndimage.generate_binary_structure(3,1) # 3D von Neumann neighborhood
for smoothIter in xrange(options.N): for smoothIter in xrange(options.N):
boundary = numpy.zeros(microstructure.shape) boundary = np.zeros(microstructure.shape)
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):
interfaceEnergy = numpy.maximum(boundary, interfaceEnergy = np.maximum(boundary,
interfacialEnergy(microstructure,numpy.roll(numpy.roll(numpy.roll( interfacialEnergy(microstructure,np.roll(np.roll(np.roll(
microstructure,i,axis=0), j,axis=1), k,axis=2))) # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood) microstructure,i,axis=0), j,axis=1), k,axis=2))) # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood)
periodic_interfaceEnergy = numpy.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,
grid[1]/2:-grid[1]/2, grid[1]/2:-grid[1]/2,
grid[2]/2:-grid[2]/2] # periodically extend interfacial energy array by half a grid size in positive and negative directions grid[2]/2:-grid[2]/2] # periodically extend interfacial energy array by half a grid size in positive and negative directions
index = ndimage.morphology.distance_transform_edt(periodic_interfaceEnergy == 0., # transform bulk volume (i.e. where interfacial energy is zero) index = ndimage.morphology.distance_transform_edt(periodic_interfaceEnergy == 0., # transform bulk volume (i.e. where interfacial energy is zero)
return_distances = False, return_distances = False,
return_indices = True) # want array index of nearest voxel on periodically extended boundary return_indices = True) # want array index of nearest voxel on periodically extended boundary
# boundaryExt = boundaryExt[index[0].flatten(),index[1].flatten(),index[2].flatten()].reshape(boundaryExt.shape) # fill bulk with energy of nearest interface | question PE: what "flatten" for? # boundaryExt = boundaryExt[index[0].flatten(),index[1].flatten(),index[2].flatten()].reshape(boundaryExt.shape) # fill bulk with energy of nearest interface | question PE: what "flatten" for?
periodic_bulkEnergy = periodic_interfaceEnergy[index[0], periodic_bulkEnergy = periodic_interfaceEnergy[index[0],
index[1], index[1],
index[2]].reshape(2*grid) # fill bulk with energy of nearest interface index[2]].reshape(2*grid) # fill bulk with energy of nearest interface
diffusedEnergy = numpy.fft.irfftn(numpy.fft.rfftn(numpy.where(ndimage.morphology.binary_dilation(interfaceEnergy > 0., diffusedEnergy = np.fft.irfftn(np.fft.rfftn(np.where(ndimage.morphology.binary_dilation(interfaceEnergy > 0.,
structure = struc, structure = struc,
iterations = options.d/2 + 1), # fat boundary | question PE: why 2d - 1? I would argue for d/2 + 1 !! iterations = options.d/2 + 1), # fat boundary | question PE: why 2d - 1? I would argue for d/2 + 1 !!
periodic_bulkEnergy[grid[0]/2:-grid[0]/2, # retain filled energy on fat boundary... periodic_bulkEnergy[grid[0]/2:-grid[0]/2, # retain filled energy on fat boundary...
@ -207,7 +189,7 @@ for file in files:
grid[2]/2:-grid[2]/2], # ...and zero everywhere else grid[2]/2:-grid[2]/2], # ...and zero everywhere else
0.)\ 0.)\
)*gauss) )*gauss)
periodic_diffusedEnergy = numpy.tile(diffusedEnergy,(3,3,3))[grid[0]/2:-grid[0]/2, periodic_diffusedEnergy = np.tile(diffusedEnergy,(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 smoothed bulk energy grid[2]/2:-grid[2]/2] # periodically extend the smoothed bulk energy
index = ndimage.morphology.distance_transform_edt(periodic_diffusedEnergy >= 0.5, # transform voxels close to interface region | question PE: what motivates 1/2 (could be any small number, or)? index = ndimage.morphology.distance_transform_edt(periodic_diffusedEnergy >= 0.5, # transform voxels close to interface region | question PE: what motivates 1/2 (could be any small number, or)?
@ -219,21 +201,21 @@ for file in files:
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
immutable = numpy.zeros(microstructure.shape, dtype=bool) immutable = np.zeros(microstructure.shape, dtype=bool)
for micro in options.immutable: for micro in options.immutable:
immutable += numpy.logical_or(microstructure == micro, microstructure_original == micro) # find locations where immutable microstructures have been or are now immutable += np.logical_or(microstructure == micro, microstructure_original == micro) # find locations where immutable microstructures have been or are now
microstructure = numpy.where(immutable, microstructure_original,microstructure) # undo any changes involving immutable microstructures microstructure = np.where(immutable, microstructure_original,microstructure) # undo any changes involving immutable microstructures
# --- renumber to sequence 1...Ngrains if requested ------------------------------------------------ # --- renumber to sequence 1...Ngrains if requested ------------------------------------------------
# http://stackoverflow.com/questions/10741346/numpy-frequency-counts-for-unique-values-in-an-array # http://stackoverflow.com/questions/10741346/np-frequency-counts-for-unique-values-in-an-array
if options.renumber: if options.renumber:
newID = 0 newID = 0
for microstructureID,count in enumerate(numpy.bincount(microstructure.flatten())): for microstructureID,count in enumerate(np.bincount(microstructure.flatten())):
if count != 0: if count != 0:
newID += 1 newID += 1
microstructure = numpy.where(microstructure == microstructureID, newID, microstructure) microstructure = np.where(microstructure == microstructureID, newID, microstructure)
# --- assemble header ----------------------------------------------------------------------------- # --- assemble header -----------------------------------------------------------------------------
newInfo['microstructures'] = microstructure[0:info['grid'][0],0:info['grid'][1],0:info['grid'][2]].max() newInfo['microstructures'] = microstructure[0:info['grid'][0],0:info['grid'][1],0:info['grid'][2]].max()
@ -257,11 +239,11 @@ for file in files:
# --- write microstructure information ------------------------------------------------------------ # --- write microstructure information ------------------------------------------------------------
formatwidth = int(math.floor(math.log10(microstructure.max())+1)) formatwidth = int(math.floor(math.log10(microstructure.max())+1))
theTable.data = microstructure[0:info['grid'][0],0:info['grid'][1],0:info['grid'][2]].reshape(numpy.prod(info['grid']),order='F').transpose() # question PE: this assumes that only the Z dimension can be 1! theTable.data = microstructure[0:info['grid'][0],0:info['grid'][1],0:info['grid'][2]].reshape(np.prod(info['grid']),order='F').transpose() # question PE: this assumes that only the Z dimension can be 1!
theTable.data_writeArray('%%%ii'%(formatwidth),delimiter=' ') theTable.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
#--- output finalization -------------------------------------------------------------------------- #--- output finalization --------------------------------------------------------------------------
if file['name'] != 'STDIN': if file['name'] != 'STDIN':
theTable.input_close() theTable.__IO__['in'].close()
theTable.output_close() theTable.__IO__['out'].close()
os.rename(file['name']+'_tmp',file['name']) os.rename(file['name']+'_tmp',file['name'])