legacy code geom... can be achieved by combining voronoi tesselation and geom_rescale, for ipfs can be plotted from ascii table with matlab

This commit is contained in:
Martin Diehl 2015-10-08 17:41:35 +00:00
parent 9ec7c08bb6
commit fd04e0763b
2 changed files with 0 additions and 403 deletions

View File

@ -1,250 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import threading,time,os,subprocess,shlex,string
import os,re,sys,math,string
import numpy as np
from optparse import OptionParser
import damask
from collections import defaultdict
scriptID = string.replace('$Id: geom_directionalNonEquiaxedGrain.py 4290 2015-07-24 08:41:08Z hm.zhang $','\n','\\n')
scriptName = os.path.splitext(scriptID.split()[1])[0]
def execute(cmd,streamIn=None,wd='./'):
'''
executes a command in given directory and returns stdout and stderr for optional stdin
'''
initialPath=os.getcwd()
os.chdir(wd)
process = subprocess.Popen(shlex.split(cmd),stdout=subprocess.PIPE,stderr = subprocess.PIPE,stdin=subprocess.PIPE)
if streamIn != None:
out,error = process.communicate(streamIn.read())
else:
out,error = process.communicate()
os.chdir(initialPath)
return out,error
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
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),
'microstructures': lambda x: int(x),
}
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Generate the geometry description of a directional non-equiaxed grain structure, e.g., RVE cutted from a cold-rolling sheet.
The initial equiaxed grain structure is generated by standard Voronoi tessellation, '--reduct' specifies the thickness
reduction after rolling, and '-n' specified the number of samples cutted, e.g., n=5, then five samples will be cutted from the
sheet along \\theta = 0 (the rolling direction), \\theta = 22.5, \\theta = 45, \\theta = 67.5, and \\theta = 90 (the
transversal direction ).
""", version = scriptID)
parser.add_option('-N', dest='N', type='int', metavar='int',
help='number of seed points to distribute [%default]')
parser.add_option('-r', '--rnd', dest='randomSeed', type='int', nargs = 2, metavar=' '.join(['int']*2),
help='seed of random number generator [%default]')
parser.add_option('-m', '--microstructure', dest='microstructure', type='int', metavar='int',
help='first microstructure index [%default]')
parser.add_option('-g', '--grid', dest='grid', type='int', nargs = 3, metavar=' '.join(['int']*3),
help='a,b,c grid of hexahedral box [from seeds file]')
parser.add_option('-s', '--size', dest='size', type='float', nargs = 3, metavar=' '.join(['float']*3),
help='x,y,z size of hexahedral box [1.0 along largest grid point number]')
parser.add_option('--phase', dest='phase', type='int', metavar = 'int',
help='phase index to be used [%default]')
parser.add_option('--crystallite', dest='crystallite', type='int', metavar = 'int',
help='crystallite index to be used [%default]')
parser.add_option('-c', '--configuration', dest='config', action='store_true',
help='output material configuration [%default]')
parser.add_option('--secondphase', type='float', dest='secondphase', metavar= 'float',
help='volume fraction of randomly distribute second phase [%default]')
parser.add_option('-l', '--laguerre', dest='laguerre', action='store_true',
help='use Laguerre (weighted Voronoi) tessellation [%default]')
parser.add_option('-n', dest='number', type='int', metavar='int',
help='the angle(degree) between the longitudinal direction of RVE and the rolling direction [%default]')
parser.add_option('--reduct', dest='reduction', type='float', metavar='float',
help='thickness reduction of rolling [%default]')
parser.set_defaults(
N = 500,
grid = (200,100,50),
size = (2.0,1.0,0.5),
phase = 1,
crystallite = 1,
secondphase = 0.0,
microstructure = 1,
laguerre = False,
randomSeed = (None,None),
config = False,
number = 5,
reduction = 0.6
)
(options,filenames) = parser.parse_args()
options.grid = np.array(options.grid)
sizeX = sizeY = max(options.size[0], options.size[1])
gridX, gridY = int(np.ceil(sizeX/options.size[0]*options.grid[0]))+1,int(np.ceil(sizeX/options.size[1]*options.grid[1]))+1
gridx, gridy, gridz = options.grid; sizex, sizey, sizez = options.size
nGrids = gridx*gridy*gridz; avgGrids = nGrids/options.N
dx, dy = options.size[0]/options.grid[0], options.size[1]/options.grid[1]
Ngrains = int(np.ceil(sizeX*sizeY/options.size[0]/options.size[1]*options.N))
filename = 'grains'+str(Ngrains)+'_'+str(gridX)+str(gridY)+str(options.grid[2])
thickness = 1.0-options.reduction
print 'run seeds_fromRandom'
execute('seeds_fromRandom -N %i -g %i %i %i %s.seeds'%(Ngrains,int(gridX*thickness)+1, gridY, int(gridz/thickness)+1, filename))
print 'run geom_fromVoronoiTessellation'
execute('geom_fromVoronoiTessellation -s %s %s %s < %s.seeds'%(sizeX*thickness, sizeY, sizez/thickness,filename))
print 'run geom_rescale'
execute('geom_rescale -g %i %i %i -s %s %s %s < %s.geom > %s_scale.geom'%(gridX, gridY, gridz,
sizeX, sizeY, sizez,filename,filename) )
print ('the size of the cutted RVE is %sX%sX%s'%(sizex, sizey, sizez) )
print ('the thickness reduction is %s'%('{:.1%}'.format(options.reduction)))
# --- loop over input files -------------------------------------------------------------------------
filenames = [filename+'.geom']
for theta in np.linspace(0, np.pi/2, options.number):
postfix = str(int(np.round(theta*180.0/np.pi)))
c_t, s_t = np.cos(theta), np.sin(theta)
offsetX, offsetY = 0.5*( sizeX - (sizex*c_t-sizey*s_t) ), 0.5*( sizeY - (sizex*s_t+sizey*c_t) )
for name in filenames:
if name == 'STDIN':
file = {'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout, 'croak':sys.stderr}
file['croak'].write('\033[1m'+scriptName+'\033[0m\n')
else:
if not os.path.exists(name): continue
file = {'name':name, 'input':open(name), 'output':open(name+postfix+'_tmp','w'), 'croak':sys.stderr}
file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n')
table = damask.ASCIItable(file['input'],file['output'],buffered=False) # make unbuffered ASCII_table
table.head_read() # read ASCII header info
#--- interpret header ----------------------------------------------------------------------------
info = {
'grid': np.zeros(3,'i'),
'size': np.array((0.0,0.0,0.0)),
'origin': np.zeros(3,'d'),
'microstructures': 0,
'homogenization': 0,
}
newInfo = {
'microstructures': 0,
}
extra_header = []
for header in table.info:
headitems = map(str.lower,header.split())
if len(headitems) == 0: continue
if headitems[0] in mappings.keys():
if headitems[0] in identifiers.keys():
for i in xrange(len(identifiers[headitems[0]])):
info[headitems[0]][i] = \
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)
newInfo['microstructures'] = info['microstructures']
if 0 not in options.grid: # user-specified grid
info['grid'] = np.array(options.grid)
for i in xrange(3):
if info['size'][i] <= 0.0: # any invalid size?
info['size'][i] = float(info['grid'][i])/max(info['grid'])
file['croak'].write('rescaling size %s...\n'%{0:'x',1:'y',2:'z'}[i])
if np.any(info['grid'] < 1):
file['croak'].write('invalid grid a b c.\n')
continue
if np.any(info['size'] <= 0.0):
file['croak'].write('invalid size x y z.\n')
continue
# read the topological data from maternal RVE file
print ( 'cut the RVE along the direction of theta = %i'%(np.round(theta*180/np.pi)) )
GrainNo = np.chararray((gridX, gridY, gridz), itemsize=6)
for i in xrange(gridY*gridz):
content = file['input'].readline().split()
for j in xrange(gridX): GrainNo[j, np.mod(i,gridY), i/gridY] = content[j]
# cut a sub-RVE cooresponding to the specified direction from the maternal RVE
subGrainNo = np.chararray((gridx, gridz, gridy), itemsize=6)
for i in xrange(gridx):
for j in xrange(gridy):
I = int(np.floor( (dx*(i+0.5)*c_t - dy*(j+0.5)*s_t + offsetX)/dx))
J = int(np.floor( (dx*(i+0.5)*s_t + dy*(j+0.5)*c_t + offsetY)/dy))
I = min(I, gridX-1); J = min(J, gridY-1)
for k in xrange(gridz): subGrainNo[i,k,j] = GrainNo[I,J,k]
subGrainNoVec = subGrainNo.reshape(nGrids)
# count the number of grains in the sub-RVE
index = defaultdict(list)
for i in xrange(nGrids): index[subGrainNoVec[i]].append(i)
ngrains = len(index)
# count the broken (scattered) grains due to the cutting, and merge them.
if ngrains > options.N*1.1:
N1, N2 = 0, nGrids
for key in index:
if len(index[key])>=0.4*avgGrids:
N1+=1; N2-=len(index[key]) # N1: valid grains; N2: number of grids needed to be re-assigned orientation
ngrid2 = min(int(0.8*avgGrids), N2/(options.N-N1)+1) # grid in each grain
newGrains = N2/ngrid2+1
a = [ngrid2]*(newGrains)
for key in index:
if len(index[key])<0.4*avgGrids:
for i in xrange(newGrains):
if a[i] >= len(index[key]):
for j in index[key]: subGrainNoVec[j] = -i-1
a[i] -= len(index[key])
break
index = defaultdict(list)
for i in xrange(nGrids): index[subGrainNoVec[i]].append(i)
ngrains = len(index)
# assign orientations
grainsNo = np.arange(ngrains); np.random.shuffle(grainsNo)
for i,key in enumerate(index):
no = str(grainsNo[i]+1)
for j in index[key]: subGrainNoVec[int(j)] = no
newInfo['microstructures'] = ngrains
if newInfo['microstructures'] == 0:
file['croak'].write('no grain info found.\n')
continue
#--- write header ---------------------------------------------------------------------------------
table.labels_clear()
table.info_clear()
table.info_append(extra_header+[
scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta %i\tb %i\tc %i"%(info['grid'][0],info['grid'][1],info['grid'][2],),
"size\tx %f\ty %f\tz %f"%(options.size[0],options.size[1],options.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']),
])
table.head_write()
# --- write microstructure information ------------------------------------------------------------
formatwidth = 1+int(math.log10(newInfo['microstructures']))
table.data = np.array(map(int, subGrainNoVec)).reshape(gridx, gridy*gridz).T
table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
#--- output finalization --------------------------------------------------------------------------
if file['name'] != 'STDIN':
prefix = os.path.splitext(file['name'])[0].replace('grains'+str(Ngrains),'grains'+str(ngrains))
os.rename(name+postfix+'_tmp',
prefix+'%s'%('_material.config' if options.config else '_'+postfix+'.geom'))
table.close()

View File

@ -1,153 +0,0 @@
#!/usr/bin/env python
import damask
import os,sys,math,re,string
from optparse import OptionParser
scriptID = string.replace('$Id$','\n','\\n')
scriptName = scriptID.split()[1]
def integerFactorization(i):
j = int(math.floor(math.sqrt(float(i))))
while (j>1 and int(i)%j != 0):
j -= 1
return j
def positiveRadians(angle):
angle = math.radians(float(angle))
while angle < 0.0:
angle += 2.0*math.pi
return angle
def getHeader(sizeX,sizeY,step,format):
if format == 'ang':
return [
'# TEM_PIXperUM 1.000000',
'# x-star 0.509548',
'# y-star 0.795272',
'# z-star 0.611799',
'# WorkingDistance 18.000000',
'#',
'# Phase 1',
'# MaterialName Al',
'# Formula Fe',
'# Info',
'# Symmetry 43',
'# LatticeConstants 2.870 2.870 2.870 90.000 90.000 90.000',
'# NumberFamilies 4',
'# hklFamilies 1 1 0 1 0.000000 1',
'# hklFamilies 2 0 0 1 0.000000 1',
'# hklFamilies 2 1 1 1 0.000000 1',
'# hklFamilies 3 1 0 1 0.000000 1',
'# Categories 0 0 0 0 0 ',
'#',
'# GRID: SquareGrid',
'# XSTEP: ' + str(step),
'# YSTEP: ' + str(step),
'# NCOLS_ODD: ' + str(sizeX),
'# NCOLS_EVEN: ' + str(sizeX),
'# NROWS: ' + str(sizeY),
'#',
'# OPERATOR: ODFsammpling',
'#',
'# SAMPLEID: ',
'#',
'# SCANID: ',
'#'
]
else: # ctf format
return [ \
'Channel Text File',
'Prj X:\\xxx\\xxxx.cpr',
'Author [MPIE-DAMASK]',
'JobMode Grid',
'XCells '+ str(sizeX),
'YCells '+ str(sizeY),
'XStep '+ str(step),
'YStep '+ str(step),
'AcqE1 0',
'AcqE2 90',
'AcqE3 0',
'Euler angles refer to Sample Coordinate system (CS0)! Mag 300 Coverage 100 Device 0 KV 20 TiltAngle 70 TiltAxis 0',
'Phases 1',
'4.05;4.05;4.05 90;90;90 Aluminium 11 225 3803863129_5.0.6.3 -2102160418 Cryogenics18,54-55',
'Phase X Y Bands Error Euler1 Euler2 Euler3 MAD BC BS'
]
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
output the Euler angles in the format supported by TSL (.ang or .ctf).
""", version = scriptID)
parser.add_option("-c", "--column", type="int", dest="column",
help="starting column of Euler triplet")
parser.add_option("-s", "--skip", type="int", dest="skip",
help="skip this many lines of heading info [%default]")
parser.add_option("-f", "--format", type="string", dest="format",
help="the format of the output file [%default]")
parser.set_defaults (column = 1)
parser.set_defaults (skip = 1)
parser.set_defaults (format = 'ang')
(options,filenames) = parser.parse_args()
options.column -= 1
#--- setup file handles ---------------------------------------------------------------------------
files = []
if filenames == []:
files.append({'name':'STDIN','input':sys.stdin,'output':sys.stdout,'croak':sys.stderr})
else:
for name in filenames:
if os.path.exists(name):
files.append({'name':name,'input':open(name),'output':open(name+'_tmp','w'),'croak':sys.stdout})
else:
print('No such file or directory: '+name)
#--- loop over input files ------------------------------------------------------------------------
for file in files:
file['croak'].write('\033[1m' + scriptName + '\033[0m: ' + (file['name'] if file['name'] != 'STDIN' else '') + '\n')
# open texture file and read content
textureFile = open(file['name'])
content = textureFile.readlines()
textureFile.close()
m = re.match('(\d+)\s+head',content[0],re.I)
if m != None and options.skip == 0: options.skip = int(m.group(1))+1
# extract orientation angles
if options.format == 'ang': # 'ang' file, radian
angles = [map(positiveRadians,line.split()[options.column:options.column+3]) for line in content[options.skip:]]
else: # 'ctf' file, degree
angles = [line.split()[options.column:options.column+3] for line in content[options.skip:]]
nPoints = len(angles)
sizeY = integerFactorization(nPoints)
sizeX = nPoints / sizeY
file['croak'].write('%s: %i*%i = %i (== %i)\n'%(file['name'],sizeX,sizeY,sizeX*sizeY,nPoints) )
# write ang/ctf file
for line in getHeader(sizeX,sizeY,1.0,options.format):
file['output'].write(line + '\n')
for counter,point in enumerate(angles):
if options.format == 'ang':
file['output'].write(''.join(['%10.5f'%angle for angle in point])+
''.join(['%10.5f'%coord for coord in [counter%sizeX,counter//sizeX]])+
' 100.0 1.0 0 1 1.0\n')
else:
file['output'].write('1 '+
' '.join(['%6.1f'%coord for coord in [counter%sizeX,counter//sizeX]])+
' 5 0 '+
' '.join(['%10.5f'%float(angle) for angle in point])+
' 0.5000 100 0\n')
if file['name'] != 'STDIN':
file['output'].close()
os.rename(file['name']+'_tmp', os.path.splitext(file['name'])[0]+'.'+options.format )