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

bug fixes and more optimizations

This commit is contained in:
Pratheek Shanthraj 2013-07-03 00:48:21 +00:00
parent 8f146ad385
commit 930e605afc
1 changed files with 51 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
processing/pre/geom_grainGrowth.py
View File

@ -27,54 +27,75 @@ class extendedOption(Option):
def grainCoarsenLocal(microLocal,ix,iy,iz,window): def grainCoarsenLocal(microLocal,ix,iy,iz,window):
interfacialEnergy = lambda A,B: 1.0 interfacialEnergy = lambda A,B: 1.0
winner = numpy.zeros(microLocal.shape).astype(int) struc = ndimage.generate_binary_structure(3,3)
winner = numpy.where(numpy.reshape(numpy.in1d(microLocal,options.black),microLocal.shape), winner = numpy.where(numpy.in1d(microLocal,options.black).reshape(microLocal.shape),
microLocal,0) # zero out non-blacklisted microstructures microLocal,0) # zero out non-blacklisted microstructures
diffusedMax = (winner > 0).astype(float) # concentration of immutable microstructures diffusedMax = (winner > 0).astype(float) # concentration of immutable microstructures
boundingSlice = ndimage.measurements.find_objects(microLocal) # bounding boxes of each distinct microstructure region boundingSlice = ndimage.measurements.find_objects(microLocal) # bounding boxes of each distinct microstructure region
for grain in set(numpy.unique(microLocal)) - set(options.black) - (set([0])): # diffuse all microstructures except immutable ones diffused = {}
for grain in set(numpy.unique(microLocal)) - set(options.black) - (set([0])): # all microstructures except immutable ones
mini = [max(0, boundingSlice[grain-1][i].start - window) for i in range(3)] # upper right of expanded bounding box mini = [max(0, boundingSlice[grain-1][i].start - window) for i in range(3)] # upper right of expanded bounding box
maxi = [min(microLocal.shape[i], boundingSlice[grain-1][i].stop + window) for i in range(3)] # lower left of expanded bounding box maxi = [min(microLocal.shape[i], boundingSlice[grain-1][i].stop + window) for i in range(3)] # lower left of expanded bounding box
microWindow = microLocal[mini[0]:maxi[0],mini[1]:maxi[1],mini[2]:maxi[2]] microWindow = microLocal[mini[0]:maxi[0],mini[1]:maxi[1],mini[2]:maxi[2]]
neighbours = set(numpy.unique(ndimage.morphology.binary_dilation(microWindow==grain,\ grainCharFunc = microWindow==grain
structure=ndimage.generate_binary_structure(3,3))))\ neighbours = set(numpy.unique(microWindow[ndimage.morphology.binary_dilation(grainCharFunc,structure=struc)]))\
- set([grain]) - set(options.black) # who is on my boundary? - set([grain]) - set(options.black) # who is on my boundary?
if len(neighbours) == 0: # no neighbours try: # has grain been diffused previously?
diffused = numpy.ones(microWindow.shape) diff = diffused[grain].copy() # yes: use previously diffused grain
elif len(neighbours) == 1: # 1 neighbour except:
speed = interfacialEnergy(grain,neighbours.pop()) # speed proportional to interfacial energy between me and only neighbour diffused[grain] = ndimage.filters.gaussian_filter((grainCharFunc).astype(float),options.d) # no: diffuse grain with unit speed
diffused = ndimage.filters.gaussian_filter((microWindow==grain).astype(float),speed*options.d)# diffuse microstructure inside extended bounding box diff = diffused[grain].copy()
else: # more than 1 neighbour ie junctions if len(neighbours) == 1:
speed = interfacialEnergy(grain,neighbours.pop()) # speed proportional to interfacial energy between me and neighbour
diff = speed*diff + (1.-speed)*(grainCharFunc) # rescale diffused microstructure by speed
else:
tiny = 0.001
numerator = numpy.zeros(microWindow.shape) numerator = numpy.zeros(microWindow.shape)
denominator = numpy.zeros(microWindow.shape) denominator = numpy.zeros(microWindow.shape)
diffused = ndimage.filters.gaussian_filter((microWindow==grain).astype(float),options.d) # diffuse microstructure inside extended bounding box weights = {grain: diff + tiny}
weights = {}
weights[grain] = diffused
for i in neighbours: for i in neighbours:
weights[i] = ndimage.filters.gaussian_filter((microWindow==i).astype(float),options.d) # partition of unity around me miniI = [max(0, boundingSlice[i-1][j].start - window) for j in range(3)] # bounding box of neighbour
for grainA,grainB in itertools.combinations(neighbours,2): # combinations of triple junctions possible maxiI = [min(microLocal.shape[j], boundingSlice[i-1][j].stop + window) for j in range(3)]
miniCommon = [max(mini[j], boundingSlice[i-1][j].start - window) for j in range(3)] # intersection of mine and neighbouring bounding box
maxiCommon = [min(maxi[j], boundingSlice[i-1][j].stop + window) for j in range(3)]
weights[i] = tiny*numpy.ones(microWindow.shape)
try: # has neighbouring grain been diffused previously?
weights[i][miniCommon[0] - mini[0]:maxiCommon[0] - mini[0],\
miniCommon[1] - mini[1]:maxiCommon[1] - mini[1],\
miniCommon[2] - mini[2]:maxiCommon[2] - mini[2]] = diffused[i][miniCommon[0] - miniI[0]:maxiCommon[0] - miniI[0],\
miniCommon[1] - miniI[1]:maxiCommon[1] - miniI[1],\
miniCommon[2] - miniI[2]:maxiCommon[2] - miniI[2]] + tiny
except: # no: diffuse neighbouring grain with unit speed
diffused[i] = ndimage.filters.gaussian_filter((microLocal[miniI[0]:maxiI[0],\
miniI[1]:maxiI[1],\
miniI[2]:maxiI[2]]==i).astype(float),options.d)
weights[i][miniCommon[0] - mini[0]:maxiCommon[0] - mini[0],\
miniCommon[1] - mini[1]:maxiCommon[1] - mini[1],\
miniCommon[2] - mini[2]:maxiCommon[2] - mini[2]] = diffused[i][miniCommon[0] - miniI[0]:maxiCommon[0] - miniI[0],\
miniCommon[1] - miniI[1]:maxiCommon[1] - miniI[1],\
miniCommon[2] - miniI[2]:maxiCommon[2] - miniI[2]] + tiny
for grainA,grainB in itertools.combinations(neighbours,2): # combinations of possible triple junctions
speed = interfacialEnergy(grain,grainA) +\ speed = interfacialEnergy(grain,grainA) +\
interfacialEnergy(grain,grainB) -\ interfacialEnergy(grain,grainB) -\
interfacialEnergy(grainA,grainB) # speed of the triple junction interfacialEnergy(grainA,grainB) # speed of the triple junction
weight = weights[grain] + weights[grainA] + weights[grainB] weight = weights[grainA] * weights[grainB]
numerator += weight*(speed*diffused + (1.-speed)*(microWindow==grain)) if numpy.any(weight > 0.01): # strongly interacting triple junction?
denominator += weight weight *= weights[grain]
numerator += weight*(speed*diff + (1.-speed)*(grainCharFunc))
denominator += weight
diffused = (numerator)/(denominator+1e-30) diff = numerator/denominator # weighted sum of strongly interacting triple junctions
isMax = diffused > diffusedMax[mini[0]:maxi[0],\
mini[1]:maxi[1],\
mini[2]:maxi[2]] # me at highest concentration?
winner[mini[0]:maxi[0],\ winner[mini[0]:maxi[0],\
mini[1]:maxi[1],\ mini[1]:maxi[1],\
mini[2]:maxi[2]][isMax] = grain # remember me ... mini[2]:maxi[2]][diff > diffusedMax[mini[0]:maxi[0],\
mini[1]:maxi[1],\
mini[2]:maxi[2]]] = grain # remember me ...
diffusedMax[mini[0]:maxi[0],\ diffusedMax[mini[0]:maxi[0],\
mini[1]:maxi[1],\ mini[1]:maxi[1],\
mini[2]:maxi[2]] = numpy.where(isMax,\ mini[2]:maxi[2]] = numpy.maximum(diff,diffusedMax[mini[0]:maxi[0],\
diffused,\ mini[1]:maxi[1],\
diffusedMax[mini[0]:maxi[0],\ mini[2]:maxi[2]]) # ... and my concentration
mini[1]:maxi[1],\
mini[2]:maxi[2]]) # ... and my concentration
return [winner[window:-window,window:-window,window:-window],ix,iy,iz] return [winner[window:-window,window:-window,window:-window],ix,iy,iz]