From e635b06270f0392f381e1f679183c813df43f583 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 4 Jun 2013 12:56:57 +0000 Subject: [PATCH] last polishing on preprocessing scripts, documentation and scripts are all up to date. added addSchmidfactors to post processing scripts --- processing/post/addSchmidfactors.py | 369 ++++++++++++++++++ processing/pre/geom_canvas.py | 4 +- processing/pre/geom_fromEuclideanDistance.py | 5 +- processing/pre/geom_fromMinimalSurface.py | 6 +- .../pre/geom_fromVoronoiTessellation.py | 2 +- processing/pre/geom_vicinityOffset.py | 3 +- .../pre/patchFromReconstructedBoundaries.py | 21 +- processing/pre/randomSeeding.py | 4 +- processing/setup/symLink_Processing.py | 1 + 9 files changed, 392 insertions(+), 23 deletions(-) create mode 100755 processing/post/addSchmidfactors.py diff --git a/processing/post/addSchmidfactors.py b/processing/post/addSchmidfactors.py new file mode 100755 index 000000000..ab9be3ec5 --- /dev/null +++ b/processing/post/addSchmidfactors.py @@ -0,0 +1,369 @@ +#!/usr/bin/python + +import os,re,sys,math +from optparse import OptionParser + +CoverA=1.587 +slipnormal_temp = [ # This is the real slip system information for hex aka titanium for now. + [0,0,0,1], + [0,0,0,1], + [0,0,0,1], + [0,1,-1,0], + [-1,0,1,0], + [1,-1,0,0], + [0,1,-1,1], + [-1,1,0,1], + [-1,0,1,1], + [0,-1,1,1], + [1,-1,0,1], + [1,0,-1,1], + [0,1,-1,1], + [0,1,-1,1], + [-1,1,0,1], + [-1,1,0,1], + [-1,0,1,1], + [-1,0,1,1], + [0,-1,1,1], + [0,-1,1,1], + [1,-1,0,1], + [1,-1,0,1], + [1,0,-1,1], + [1,0,-1,1], + ] + +slipdirection_temp = [ + [2,-1,-1,0], + [-1,2,-1,0], + [-1,-1,2,0], + [2,-1,-1,0], + [-1,2,-1,0], + [-1,-1,2,0], + [2,-1,-1,0], + [1,1,-2,0], + [-1,2,-1,0], + [-2,1,1,0], + [-1,-1,2,0], + [1,-2,1,0], + [-1,2,-1,3], + [1,1,-2,3], + [-2,1,1,3], + [-1,2,-1,3], + [-1,-1,2,3], + [-2,1,1,3], + [1,-2,1,3], + [-1,-1,2,3], + [2,-1,-1,3], + [1,-2,1,3], + [1,1,-2,3], + [2,-1,-1,3], + ] + +# slip normals and directions according to cpfem implementation +Nslipsystems = {'fcc': 12, 'bcc': 24, 'hex': 24} +slipnormal = { \ + 'fcc': [ + [1,1,1], + [1,1,1], + [1,1,1], + [-1,-1,1], + [-1,-1,1], + [-1,-1,1], + [1,-1,-1], + [1,-1,-1], + [1,-1,-1], + [-1,1,-1], + [-1,1,-1], + [-1,1,-1], + ], + 'bcc': [ + [0,1,1], + [0,1,1], + [0,-1,1], + [0,-1,1], + [1,0,1], + [1,0,1], + [-1,0,1], + [-1,0,1], + [1,1,0], + [1,1,0], + [-1,1,0], + [-1,1,0], + [2,1,1], + [-2,1,1], + [2,-1,1], + [2,1,-1], + [1,2,1], + [-1,2,1], + [1,-2,1], + [1,2,-1], + [1,1,2], + [-1,1,2], + [1,-1,2], + [1,1,-2], + ], + 'hex': [ # these are dummy numbers and are recalculated based on the above hex real slip systems. + [1,1,0], + [1,1,0], + [1,0,1], + [1,0,1], + [0,1,1], + [0,1,1], + [1,-1,0], + [1,-1,0], + [-1,0,1], + [-1,0,1], + [0,-1,1], + [0,-1,1], + [2,-1,1], + [1,-2,-1], + [1,1,2], + [2,1,1], + [1,2,-1], + [1,-1,2], + [2,1,-1], + [1,2,1], + [1,-1,-2], + [2,-1,-1], + [1,-2,1], + [1,1,-2], + ], + } +slipdirection = { \ + 'fcc': [ + [0,1,-1], + [-1,0,1], + [1,-1,0], + [0,-1,-1], + [1,0,1], + [-1,1,0], + [0,-1,1], + [-1,0,-1], + [1,1,0], + [0,1,1], + [1,0,-1], + [-1,-1,0], + ], + 'bcc': [ + [1,-1,1], + [-1,-1,1], + [1,1,1], + [-1,1,1], + [-1,1,1], + [-1,-1,1], + [1,1,1], + [1,-1,1], + [-1,1,1], + [-1,1,-1], + [1,1,1], + [1,1,-1], + [-1,1,1], + [1,1,1], + [1,1,-1], + [1,-1,1], + [1,-1,1], + [1,1,-1], + [1,1,1], + [-1,1,1], + [1,1,-1], + [1,-1,1], + [-1,1,1], + [1,1,1], + ], + 'hex': [ # these are dummy numbers and are recalculated based on the above hex real slip systems. + [-1,1,1], + [1,-1,1], + [-1,-1,1], + [-1,1,1], + [-1,-1,1], + [1,-1,1], + [1,1,1], + [-1,-1,1], + [1,-1,1], + [1,1,1], + [1,1,1], + [-1,1,1], + [1,1,-1], + [1,1,-1], + [1,1,-1], + [1,-1,-1], + [1,-1,-1], + [1,-1,-1], + [1,-1,1], + [1,-1,1], + [1,-1,1], + [1,1,1], + [1,1,1], + [1,1,1], + ], + } + + +# -------------------------------------------------------------------- + +def applyEulers(phi1,Phi,phi2,x): + """ transform x given in crystal coordinates to xbar returned in lab coordinates for Euler angles phi1,Phi,phi2 """ + + eulerRot = [[ math.cos(phi1)*math.cos(phi2) - math.cos(Phi)*math.sin(phi1)*math.sin(phi2), - math.cos(phi1)*math.sin(phi2) - math.cos(Phi)*math.cos(phi2)*math.sin(phi1), math.sin(Phi)*math.sin(phi1)], \ + [ math.cos(phi2)*math.sin(phi1) + math.cos(Phi)*math.cos(phi1)*math.sin(phi2), math.cos(Phi)*math.cos(phi1)*math.cos(phi2) - math.sin(phi1)*math.sin(phi2), -math.sin(Phi)*math.cos(phi1)], \ + [ math.sin(Phi)*math.sin(phi2), math.sin(Phi)*math.cos(phi2), math.cos(Phi)]] + + xbar = [0,0,0] + if len(x) == 3: + for i in range(3): + xbar[i] = sum([eulerRot[i][j]*x[j] for j in range(3)]) + return xbar + +# -------------------------------------------------------------------- + +def normalize(x): + + norm = math.sqrt(sum([x[i]*x[i] for i in range(len(x))])) + + return [x[i]/norm for i in range(len(x))] + +# -------------------------------------------------------------------- + +def crossproduct(x,y): + + return [ + x[1]*y[2]-y[1]*x[2], + x[2]*y[0]-y[2]*x[0], + x[0]*y[1]-y[0]*x[1], + ] + +# -------------------------------------------------------------------- + + + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- +parser = OptionParser(usage='%prog [options] [file]', description = """ +Add columns listing Schmid factors (and optional trace vector of selected system) for given Euler angles. +Column headings need to have names 'phi1', 'Phi', 'phi2'. + +$Id$ +""") + +parser.add_option('-l','--lattice', dest='lattice', choices=('fcc','bcc','hex'), \ + help='key for lattice type [%default]') +parser.add_option('-d','--forcedirection', dest='forcedirection', type='int', nargs=3, \ + help='force direction in lab coordinates [%default]') +parser.add_option('-n','--stressnormal', dest='stressnormal', type='int', nargs=3, \ + help='stress plane normal in lab coordinates [%default]') +parser.add_option('-t','--trace', dest='traceplane', type='int', nargs=3, \ + help="normal (in lab coordinates) of plane on which the plane trace of the Schmid factor(s) is reported [%default]") +parser.add_option('-r','--rank', dest='rank', type='int', \ + help="report trace of r'th highest Schmid factor [%default]") +parser.set_defaults(lattice = 'fcc') +parser.set_defaults(forcedirection = [0, 0, 1]) +parser.set_defaults(stressnormal = None) +parser.set_defaults(traceplane = None) +parser.set_defaults(rank = 0) + +(options,filename) = parser.parse_args() + +options.forcedirection = normalize(options.forcedirection) +if options.stressnormal: + if abs(sum([options.forcedirection[i] * options.stressnormal[i] for i in range(3)])) < 1e-3: + options.stressnormal = normalize(options.stressnormal) + else: + parser.error('stress plane normal not orthogonal to force direction') +else: + options.stressnormal = options.forcedirection +if options.traceplane: + options.traceplane = normalize(options.traceplane) +options.rank = min(options.rank,Nslipsystems[options.lattice]) + +# read from standard input unless input file specified +if filename == []: + file = sys.stdin +elif os.path.exists(filename[0]): + file = open(filename[0]) + +# read data +content = file.readlines() +file.close() + +# get labels by either read the first row, or - if keyword header is present - the last line of the header +headerlines = 1 +m = re.search('(\d+)\s*head', content[0].lower()) +if m: + headerlines = int(m.group(1))+1 +labels = content[headerlines-1].split() +data = content[headerlines:] + +# Convert 4 Miller indices notation of hex to orthogonal 3 Miller indices notation +if options.lattice=="hex": + for i in range(Nslipsystems[options.lattice]): + slipnormal[options.lattice][i][0]=slipnormal_temp[i][0] + slipnormal[options.lattice][i][1]=(slipnormal_temp[i][0]+2.0*slipnormal_temp[i][1])/math.sqrt(3.0) + slipnormal[options.lattice][i][2]=slipnormal_temp[i][3]/CoverA + slipdirection[options.lattice][i][0]=slipdirection_temp[i][0]*1.5 # direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] , + slipdirection[options.lattice][i][1]=(slipdirection_temp[i][0]+2.0*slipdirection_temp[i][1])*(0.5*math.sqrt(3.0)) + slipdirection[options.lattice][i][2]=slipdirection_temp[i][3]*CoverA + + for i in range(Nslipsystems[options.lattice]): + slipnormal[options.lattice][i]=normalize(slipnormal[options.lattice][i]) + slipdirection[options.lattice][i]=normalize(slipdirection[options.lattice][i]) + +for c in range(len(labels)): + m = re.search('.*([Pp]hi\d*).*', labels[c]) + if m: + if m.group(1).lower() == "phi1": + phi1Column = c + elif m.group(1).lower() == "phi": + PhiColumn = c + elif m.group(1).lower() == "phi2": + phi2Column = c + +output = '1\theader\n' + \ + '\t'.join(map(str,labels)) + \ + '\t' + \ + '\t'.join(['(%i)S(%i %i %i)[%i %i %i]'%(i+1, + slipnormal[options.lattice][i][0], + slipnormal[options.lattice][i][1], + slipnormal[options.lattice][i][2], + slipdirection[options.lattice][i][0], + slipdirection[options.lattice][i][1], + slipdirection[options.lattice][i][2], + ) for i in range(Nslipsystems[options.lattice])]) +if options.traceplane: + if options.rank > 0: + output += '\ttrace_x\ttrace_y\ttrace_z\tsystem' + else: + output += '\t' + '\t'.join(['(%i)tx\tty\ttz'%(i+1) for i in range(Nslipsystems[options.lattice])]) +output += '\n' + +for line in data: + items = line.split()[:len(labels)] + if items == []: + continue + phi1 = math.radians(float(items[phi1Column])) + Phi = math.radians(float(items[PhiColumn])) + phi2 = math.radians(float(items[phi2Column])) + S = [ sum( [applyEulers(phi1,Phi,phi2,normalize( slipnormal[options.lattice][slipsystem]))[i]*options.stressnormal[i] for i in range(3)] ) * \ + sum( [applyEulers(phi1,Phi,phi2,normalize(slipdirection[options.lattice][slipsystem]))[i]*options.forcedirection[i] for i in range(3)] ) \ + for slipsystem in range(Nslipsystems[options.lattice]) ] + output += '\t'.join(items + map(str,S)) + if options.traceplane: + trace = [crossproduct(options.traceplane,applyEulers(phi1,Phi,phi2,normalize(slipnormal[options.lattice][slipsystem]))) \ + for slipsystem in range(Nslipsystems[options.lattice]) ] + if options.rank == 0: + output += '\t' + '\t'.join(map(lambda x:'%f\t%f\t%f'%(x[0],x[1],x[2]),trace)) + elif options.rank > 0: + SabsSorted = sorted([(abs(S[i]),i) for i in range(len(S))]) + output += '\t' + '\t'.join(map(str,trace[SabsSorted[-options.rank][1]])) + '\t%i'%(1+SabsSorted[-options.rank][1]) + # for t in [normalize(crossproduct(options.traceplane,applyEulers(phi1,Phi,phi2,normalize(slipnormal[options.lattice][i])))) for i in range(12,24)]: + # print '\t'.join(map(str,t)) + # print '\t'.join(map(lambda x: str(-x),t)) + # print '\t'.join(['0','0','0']) + # print + output += '\n' + +if filename == []: + print output +else: + file = open(filename[0],'w') + file.write(output) + file.close() diff --git a/processing/pre/geom_canvas.py b/processing/pre/geom_canvas.py index 64f4741f2..fac1ac0f4 100755 --- a/processing/pre/geom_canvas.py +++ b/processing/pre/geom_canvas.py @@ -46,9 +46,9 @@ Changes the (three-dimensional) canvas of a spectral geometry description. parser.add_option('-g', '--grid', dest='grid', type='int', nargs = 3, \ help='a,b,c grid of hexahedral box [unchanged]') parser.add_option('-o', '--offset', dest='offset', type='int', nargs = 3, \ - help='x,y,z offset from old to new origin of grid %default') + help='a,b,c offset from old to new origin of grid %default') parser.add_option('-f', '--fill', dest='fill', type='int', \ - help='(background) canvas grain index [autodetect]') + help='(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]') parser.add_option('-2', '--twodimensional', dest='twoD', action='store_true', \ help='output geom file with two-dimensional data arrangement [%default]') diff --git a/processing/pre/geom_fromEuclideanDistance.py b/processing/pre/geom_fromEuclideanDistance.py index 3e5726b87..7211076bb 100755 --- a/processing/pre/geom_fromEuclideanDistance.py +++ b/processing/pre/geom_fromEuclideanDistance.py @@ -117,7 +117,7 @@ boundaries, triple lines, and quadruple points. parser.add_option('-t','--type', dest='type', action='extend', type='string', \ help='feature type (%s)'%(', '.join(map(lambda x:', '.join(x['names']),features)))) -parser.add_option('-n','--neighborhood', dest='neigborhood', action='store', type='string', \ +parser.add_option('-n','--neighborhood', dest='neigborhood', choices=neighborhoods.keys(), \ help='type of neighborhood (%s) [neumann]'%(', '.join(neighborhoods.keys()))) parser.add_option('-2', '--twodimensional', dest='twoD', action='store_true', \ help='output geom file with two-dimensional data arrangement [%default]') @@ -128,9 +128,6 @@ parser.set_defaults(twoD = False) (options,filenames) = parser.parse_args() -options.neighborhood = options.neighborhood.lower() -if options.neighborhood not in neighborhoods: - parser.error('unknown neighborhood %s!'%options.neighborhood) feature_list = [] for i,feature in enumerate(features): diff --git a/processing/pre/geom_fromMinimalSurface.py b/processing/pre/geom_fromMinimalSurface.py index cccddd438..7068c0c23 100755 --- a/processing/pre/geom_fromMinimalSurface.py +++ b/processing/pre/geom_fromMinimalSurface.py @@ -40,8 +40,8 @@ Generate a geometry file of a bicontinuous structure of given type. """ + string.replace('$Id$','\n','\\n') ) -parser.add_option('-t','--type', dest='type', type='string', \ - help='type of minimal surface (%s)'%(','.join(minimal_surfaces))) +parser.add_option('-t','--type', dest='type', choices=minimal_surfaces, \ + help='type of minimal surface (%s) [primitive]' %(','.join(minimal_surfaces))) parser.add_option('-f','--threshold', dest='threshold', type='float', \ help='threshold value defining minimal surface [%default]') parser.add_option('-g', '--grid', dest='grid', type='int', nargs=3, \ @@ -51,7 +51,7 @@ parser.add_option('-s', '--size', dest='size', type='float', nargs=3, \ parser.add_option('-p', '--periods', dest='periods', type='int', \ help='number of repetitions of unit cell [%default]') parser.add_option('--homogenization', dest='homogenization', type='int', \ - help='homogenization index to be used [%defaults]') + help='homogenization index to be used [%default]') parser.add_option('--m', dest='microstructure', type='int', nargs = 2, \ help='two microstructure indices to be used %default') parser.add_option('-2', '--twodimensional', dest='twoD', action='store_true', \ diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index d6b03915d..99f573eaf 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -40,7 +40,7 @@ Generate geometry description and material configuration by standard Voronoi tes ) parser.add_option('-g', '--grid', dest='grid', type='int', nargs = 3, \ - help='a,b,c grid of hexahedral box [from seed file]') + help='a,b,c grid of hexahedral box [from seeds file]') parser.add_option('-s', '--size', dest='size', type='float', nargs = 3, \ help='x,y,z size of hexahedral box [1.0 along largest grid point number]') parser.add_option('--homogenization', dest='homogenization', type='int', \ diff --git a/processing/pre/geom_vicinityOffset.py b/processing/pre/geom_vicinityOffset.py index 9d7f9cfca..5072f808f 100755 --- a/processing/pre/geom_vicinityOffset.py +++ b/processing/pre/geom_vicinityOffset.py @@ -48,7 +48,8 @@ i.e. within the region close to a grain/phase boundary. parser.add_option('-v', '--vicinity', dest='vicinity', type='int', \ help='voxel distance checked for presence of other microstructure [%default]') parser.add_option('-m', '--microstructureoffset', dest='offset', type='int', \ - help='integer offset for tagged microstructure [autodetect]') + help='offset (positive or negative) for tagged microstructure. '+ + '"0" selects maximum microstructure index [%default]') parser.add_option('-2', '--twodimensional', dest='twoD', action='store_true', \ help='output geom file with two-dimensional data arrangement') diff --git a/processing/pre/patchFromReconstructedBoundaries.py b/processing/pre/patchFromReconstructedBoundaries.py index 2ff42f4c4..e5e788e1c 100755 --- a/processing/pre/patchFromReconstructedBoundaries.py +++ b/processing/pre/patchFromReconstructedBoundaries.py @@ -789,10 +789,10 @@ reconstructed boundary file ) parser.add_option("-o", "--output", action='extend', dest='output', type='string', \ - help="types of output [image,mentat,procedure,spectral]") + help="types of output [image, mentat, procedure, spectral]") parser.add_option("-p", "--port", type="int",\ dest="port",\ - help="Mentat connection port") + help="Mentat connection port [%default]") parser.add_option("-2", "--twodimensional", action="store_true",\ dest="twoD",\ help="twodimensional model [%default]") @@ -804,16 +804,16 @@ parser.add_option("-e", "--strain", type="float",\ help="final strain to reach in simulation [%default]") parser.add_option("--rate", type="float",\ dest="strainrate",\ - help="(engineering) strain rate to simulate") + help="(engineering) strain rate to simulate [%default]") parser.add_option("-N", "--increments", type="int",\ dest="increments",\ - help="number of increments to take") + help="number of increments to take [%default]") parser.add_option("-t", "--tolerance", type="float",\ dest="tolerance",\ - help="relative tolerance of pixel positions to be swept") + help="relative tolerance of pixel positions to be swept [%default]") parser.add_option("-m", "--mesh", choices=['dt_planar_trimesh','af_planar_trimesh','af_planar_quadmesh'],\ dest="mesh",\ - help="algorithm and element type for automeshing [%default]") + help="algorithm and element type for automeshing (dt_planar_trimesh, af_planar_trimesh, af_planar_quadmesh) [%default]") parser.add_option("-x", "--xmargin", type="float",\ dest="xmargin",\ help="margin in x in units of patch size [%default]") @@ -828,22 +828,23 @@ parser.add_option("-z", "--extrusion", type="int",\ help="number of repetitions in z-direction [%default]") parser.add_option("-i", "--imagesize", type="int",\ dest="imgsize",\ - help="size of PNG image") + help="size of PNG image [%default]") parser.add_option("-M", "--coordtransformation", type="float", nargs=4, \ dest="M",\ - help="2x2 transformation from rcb to Euler coords ( = M . [x_rcb,y_rcb])") + help="2x2 transformation from rcb to Euler coords [%default]") parser.add_option("--scatter", type="float",\ dest="scatter",\ - help="orientation scatter [%default]") + help="orientation scatter %default") parser.set_defaults(output = []) parser.set_defaults(size = 1.0) +parser.set_defaults(port = 40007) parser.set_defaults(xmargin = 0.0) parser.set_defaults(ymargin = 0.0) parser.set_defaults(resolution = 64) parser.set_defaults(extrusion = 2) parser.set_defaults(imgsize = 512) -parser.set_defaults(M = [0.0,1.0,1.0,0]) # M_11, M_12, M_21, M_22. x,y in RCB is y,x of Eulers!! +parser.set_defaults(M = [0.0,1.0,1.0,0.0]) # M_11, M_12, M_21, M_22. x,y in RCB is y,x of Eulers!! parser.set_defaults(tolerance = 1.0e-3) parser.set_defaults(scatter = 0.0) parser.set_defaults(strain = 0.2) diff --git a/processing/pre/randomSeeding.py b/processing/pre/randomSeeding.py index 37ac89cf8..ed0eea42d 100755 --- a/processing/pre/randomSeeding.py +++ b/processing/pre/randomSeeding.py @@ -35,7 +35,7 @@ mappings = { parser = OptionParser(option_class=extendedOption, usage='%prog [options]', description = """ -Distribute given number of points randomly within the three-dimensional cube [0,0,0]--[1,1,1]. +Distribute given number of points randomly within the three-dimensional cube [0.0,0.0,0.0]--[1.0,1.0,1.0]. Reports positions with random crystal orientations in seeds file format to STDOUT. """ + string.replace('$Id$','\n','\\n') ) @@ -55,7 +55,7 @@ parser.set_defaults(N = 20) Npoints = reduce(lambda x, y: x * y, options.grid) if 0 in options.grid: - file['croak'].write('invalid grid a b c.\n') + file['croak'].write('invalid grid a b c.\n') sys.exit() if options.N > Npoints: sys.stderr.write('Warning: more seeds than grid points at minimum resolution.\n') diff --git a/processing/setup/symLink_Processing.py b/processing/setup/symLink_Processing.py index f15c0b407..b108cd098 100755 --- a/processing/setup/symLink_Processing.py +++ b/processing/setup/symLink_Processing.py @@ -47,6 +47,7 @@ bin_link = { \ 'addMises.py', 'addNorm.py', 'addPK2.py', + 'addSchmidfactors.py', 'addSpectralDecomposition.py', 'addStrainTensors.py', 'averageDown.py',