2018-11-17 12:42:12 +05:30
|
|
|
#!/usr/bin/env python3
|
2015-03-05 15:35:00 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
import os
|
|
|
|
import sys
|
|
|
|
import math
|
|
|
|
import random
|
2015-03-06 00:31:31 +05:30
|
|
|
from optparse import OptionParser
|
2020-03-18 18:17:09 +05:30
|
|
|
|
2015-04-15 20:07:46 +05:30
|
|
|
import numpy as np
|
2015-03-06 00:31:31 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
import damask
|
|
|
|
|
|
|
|
|
2016-01-27 22:36:00 +05:30
|
|
|
scriptName = os.path.splitext(os.path.basename(__file__))[0]
|
|
|
|
scriptID = ' '.join([scriptName,damask.version])
|
2015-10-13 22:32:07 +05:30
|
|
|
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
# --- helper functions ---
|
2015-06-11 22:20:29 +05:30
|
|
|
def integerFactorization(i):
|
|
|
|
|
|
|
|
j = int(math.floor(math.sqrt(float(i))))
|
2015-08-08 00:33:26 +05:30
|
|
|
while j>1 and int(i)%j != 0:
|
2015-06-11 22:20:29 +05:30
|
|
|
j -= 1
|
|
|
|
return j
|
|
|
|
|
2015-03-05 15:35:00 +05:30
|
|
|
def binAsBins(bin,intervals):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""Explode compound bin into 3D bins list."""
|
2015-03-05 15:35:00 +05:30
|
|
|
bins = [0]*3
|
|
|
|
bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0]
|
|
|
|
bins[1] = (bin//intervals[2]) % intervals[1]
|
|
|
|
bins[2] = bin % intervals[2]
|
|
|
|
return bins
|
|
|
|
|
|
|
|
def binsAsBin(bins,intervals):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""Implode 3D bins into compound bin."""
|
2015-03-05 15:35:00 +05:30
|
|
|
return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2]
|
|
|
|
|
|
|
|
def EulersAsBins(Eulers,intervals,deltas,center):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""Return list of Eulers translated into 3D bins list."""
|
2015-04-15 20:07:46 +05:30
|
|
|
return [int((euler+(0.5-center)*delta)//delta)%interval \
|
|
|
|
for euler,delta,interval in zip(Eulers,deltas,intervals) \
|
|
|
|
]
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
def binAsEulers(bin,intervals,deltas,center):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""Compound bin number translated into list of Eulers."""
|
2015-03-05 15:35:00 +05:30
|
|
|
Eulers = [0.0]*3
|
|
|
|
Eulers[2] = (bin%intervals[2] + center)*deltas[2]
|
|
|
|
Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1]
|
|
|
|
Eulers[0] = (bin//(intervals[2]*intervals[1]) + center)*deltas[0]
|
|
|
|
return Eulers
|
|
|
|
|
|
|
|
def directInvRepetitions(probability,scale):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""Calculate number of samples drawn by direct inversion."""
|
2015-03-05 15:35:00 +05:30
|
|
|
nDirectInv = 0
|
2015-08-08 00:33:26 +05:30
|
|
|
for bin in range(len(probability)): # loop over bins
|
|
|
|
nDirectInv += int(round(probability[bin]*scale)) # calc repetition
|
2015-03-05 15:35:00 +05:30
|
|
|
return nDirectInv
|
|
|
|
|
|
|
|
|
2015-04-15 20:07:46 +05:30
|
|
|
# ---------------------- sampling methods -----------------------------------------------------------------------
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
# ----- efficient algorithm ---------
|
|
|
|
def directInversion (ODF,nSamples):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
|
2015-08-08 00:33:26 +05:30
|
|
|
nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
nInvSamples = 0
|
2015-04-15 20:07:46 +05:30
|
|
|
repetition = [None]*ODF['nBins']
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
scaleLower = 0.0
|
|
|
|
nInvSamplesLower = 0
|
|
|
|
scaleUpper = float(nOptSamples)
|
|
|
|
incFactor = 1.0
|
|
|
|
nIter = 0
|
|
|
|
nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper)
|
|
|
|
while (\
|
|
|
|
(scaleUpper-scaleLower > scaleUpper*1e-15 or nInvSamplesUpper < nOptSamples) and \
|
|
|
|
nInvSamplesUpper != nOptSamples \
|
2015-08-08 00:33:26 +05:30
|
|
|
): # closer match required?
|
2015-03-05 15:35:00 +05:30
|
|
|
if nInvSamplesUpper < nOptSamples:
|
|
|
|
scaleLower,scaleUpper = scaleUpper,scaleUpper+incFactor*(scaleUpper-scaleLower)/2.0
|
|
|
|
incFactor *= 2.0
|
|
|
|
nInvSamplesLower,nInvSamplesUpper = nInvSamplesUpper,directInvRepetitions(ODF['dV_V'],scaleUpper)
|
|
|
|
else:
|
|
|
|
scaleUpper = (scaleLower+scaleUpper)/2.0
|
|
|
|
incFactor = 1.0
|
|
|
|
nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper)
|
|
|
|
nIter += 1
|
2015-10-13 22:32:07 +05:30
|
|
|
damask.util.croak('%i:(%12.11f,%12.11f) %i <= %i <= %i'%(nIter,scaleLower,scaleUpper,
|
2015-08-08 00:33:26 +05:30
|
|
|
nInvSamplesLower,nOptSamples,nInvSamplesUpper))
|
2015-03-05 15:35:00 +05:30
|
|
|
nInvSamples = nInvSamplesUpper
|
|
|
|
scale = scaleUpper
|
2015-10-13 22:32:07 +05:30
|
|
|
damask.util.croak('created set of %i samples (%12.11f) with scaling %12.11f delivering %i'%(nInvSamples,
|
2015-08-08 00:33:26 +05:30
|
|
|
float(nInvSamples)/nOptSamples-1.0,
|
|
|
|
scale,nSamples))
|
|
|
|
repetition = [None]*ODF['nBins'] # preallocate and clear
|
2015-03-05 15:35:00 +05:30
|
|
|
|
2015-08-08 00:33:26 +05:30
|
|
|
for bin in range(ODF['nBins']): # loop over bins
|
|
|
|
repetition[bin] = int(round(ODF['dV_V'][bin]*scale)) # calc repetition
|
2015-03-05 15:35:00 +05:30
|
|
|
|
2016-03-02 17:13:09 +05:30
|
|
|
# build set
|
2015-03-05 15:35:00 +05:30
|
|
|
set = [None]*nInvSamples
|
|
|
|
i = 0
|
2015-04-15 20:07:46 +05:30
|
|
|
for bin in range(ODF['nBins']):
|
2015-08-08 00:33:26 +05:30
|
|
|
set[i:i+repetition[bin]] = [bin]*repetition[bin] # fill set with bin, i.e. orientation
|
|
|
|
i += repetition[bin] # advance set counter
|
2015-03-05 15:35:00 +05:30
|
|
|
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations = np.zeros((nSamples,3),'f')
|
|
|
|
reconstructedODF = np.zeros(ODF['nBins'],'f')
|
2015-03-05 15:35:00 +05:30
|
|
|
unitInc = 1.0/nSamples
|
|
|
|
for j in range(nSamples):
|
|
|
|
if (j == nInvSamples-1): ex = j
|
|
|
|
else: ex = int(round(random.uniform(j+0.5,nInvSamples-0.5)))
|
|
|
|
bin = set[ex]
|
2015-04-15 20:07:46 +05:30
|
|
|
Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations[j] = np.degrees(Eulers)
|
2015-03-05 15:35:00 +05:30
|
|
|
reconstructedODF[bin] += unitInc
|
2015-08-08 00:33:26 +05:30
|
|
|
set[ex] = set[j] # exchange orientations
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
return orientations, reconstructedODF
|
|
|
|
|
|
|
|
|
|
|
|
# ----- trial and error algorithms ---------
|
|
|
|
|
|
|
|
def MonteCarloEulers (ODF,nSamples):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
|
2015-03-05 15:35:00 +05:30
|
|
|
countMC = 0
|
|
|
|
maxdV_V = max(ODF['dV_V'])
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations = np.zeros((nSamples,3),'f')
|
|
|
|
reconstructedODF = np.zeros(ODF['nBins'],'f')
|
2015-03-05 15:35:00 +05:30
|
|
|
unitInc = 1.0/nSamples
|
|
|
|
|
|
|
|
for j in range(nSamples):
|
|
|
|
MC = maxdV_V*2.0
|
|
|
|
bin = 0
|
|
|
|
while MC > ODF['dV_V'][bin]:
|
|
|
|
countMC += 1
|
|
|
|
MC = maxdV_V*random.random()
|
2015-04-15 20:07:46 +05:30
|
|
|
Eulers = [limit*random.random() for limit in ODF['limit']]
|
|
|
|
bins = EulersAsBins(Eulers,ODF['interval'],ODF['delta'],ODF['center'])
|
|
|
|
bin = binsAsBin(bins,ODF['interval'])
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations[j] = np.degrees(Eulers)
|
2015-03-05 15:35:00 +05:30
|
|
|
reconstructedODF[bin] += unitInc
|
|
|
|
|
|
|
|
return orientations, reconstructedODF, countMC
|
|
|
|
|
|
|
|
|
|
|
|
def MonteCarloBins (ODF,nSamples):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
|
2015-03-05 15:35:00 +05:30
|
|
|
countMC = 0
|
|
|
|
maxdV_V = max(ODF['dV_V'])
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations = np.zeros((nSamples,3),'f')
|
|
|
|
reconstructedODF = np.zeros(ODF['nBins'],'f')
|
2015-03-05 15:35:00 +05:30
|
|
|
unitInc = 1.0/nSamples
|
|
|
|
|
|
|
|
for j in range(nSamples):
|
|
|
|
MC = maxdV_V*2.0
|
|
|
|
bin = 0
|
|
|
|
while MC > ODF['dV_V'][bin]:
|
|
|
|
countMC += 1
|
|
|
|
MC = maxdV_V*random.random()
|
2015-04-15 20:07:46 +05:30
|
|
|
bin = int(ODF['nBins'] * random.random())
|
|
|
|
Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations[j] = np.degrees(Eulers)
|
2015-03-05 15:35:00 +05:30
|
|
|
reconstructedODF[bin] += unitInc
|
|
|
|
|
|
|
|
return orientations, reconstructedODF
|
|
|
|
|
|
|
|
|
|
|
|
def TothVanHoutteSTAT (ODF,nSamples):
|
2019-12-10 20:50:50 +05:30
|
|
|
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations = np.zeros((nSamples,3),'f')
|
|
|
|
reconstructedODF = np.zeros(ODF['nBins'],'f')
|
2015-03-05 15:35:00 +05:30
|
|
|
unitInc = 1.0/nSamples
|
|
|
|
|
|
|
|
selectors = [random.random() for i in range(nSamples)]
|
|
|
|
selectors.sort()
|
|
|
|
indexSelector = 0
|
|
|
|
|
|
|
|
cumdV_V = 0.0
|
|
|
|
countSamples = 0
|
|
|
|
|
2015-04-15 20:07:46 +05:30
|
|
|
for bin in range(ODF['nBins']) :
|
2015-03-05 15:35:00 +05:30
|
|
|
cumdV_V += ODF['dV_V'][bin]
|
|
|
|
while indexSelector < nSamples and selectors[indexSelector] < cumdV_V:
|
2015-04-15 20:07:46 +05:30
|
|
|
Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
|
2015-08-08 00:33:26 +05:30
|
|
|
orientations[countSamples] = np.degrees(Eulers)
|
2015-03-05 15:35:00 +05:30
|
|
|
reconstructedODF[bin] += unitInc
|
|
|
|
countSamples += 1
|
|
|
|
indexSelector += 1
|
2015-08-08 00:33:26 +05:30
|
|
|
|
2015-10-13 22:32:07 +05:30
|
|
|
damask.util.croak('created set of %i when asked to deliver %i'%(countSamples,nSamples))
|
2015-03-05 15:35:00 +05:30
|
|
|
|
|
|
|
return orientations, reconstructedODF
|
|
|
|
|
|
|
|
|
2015-03-06 00:31:31 +05:30
|
|
|
# --------------------------------------------------------------------
|
|
|
|
# MAIN
|
|
|
|
# --------------------------------------------------------------------
|
2016-03-02 17:13:09 +05:30
|
|
|
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description ="""
|
2015-08-08 00:33:26 +05:30
|
|
|
Transform linear binned ODF data into given number of orientations.
|
2015-11-16 16:22:56 +05:30
|
|
|
IA: integral approximation, STAT: Van Houtte, MC: Monte Carlo
|
2015-04-15 20:07:46 +05:30
|
|
|
|
|
|
|
""", version = scriptID)
|
|
|
|
|
2015-11-16 16:22:56 +05:30
|
|
|
algorithms = ['IA', 'STAT','MC']
|
2015-08-08 00:33:26 +05:30
|
|
|
parser.add_option('-n', '--nsamples',
|
|
|
|
dest = 'number',
|
|
|
|
type = 'int', metavar = 'int',
|
|
|
|
help = 'number of orientations to be generated [%default]')
|
|
|
|
parser.add_option('-a','--algorithm',
|
|
|
|
dest = 'algorithm',
|
2015-11-16 16:22:56 +05:30
|
|
|
choices = algorithms, metavar = 'string',
|
|
|
|
help = 'sampling algorithm {%s} [IA]'%(', '.join(algorithms)))
|
2015-08-08 00:33:26 +05:30
|
|
|
parser.add_option('-p','--phase',
|
|
|
|
dest = 'phase',
|
|
|
|
type = 'int', metavar = 'int',
|
|
|
|
help = 'phase index to be used [%default]')
|
|
|
|
parser.add_option('-r', '--rnd',
|
|
|
|
dest = 'randomSeed',
|
|
|
|
type = 'int', metavar = 'int', \
|
|
|
|
help = 'seed of random number generator [%default]')
|
|
|
|
parser.set_defaults(randomSeed = None,
|
|
|
|
number = 500,
|
|
|
|
algorithm = 'IA',
|
|
|
|
phase = 1,
|
|
|
|
ang = True,
|
|
|
|
)
|
|
|
|
|
2015-04-15 20:07:46 +05:30
|
|
|
(options,filenames) = parser.parse_args()
|
2015-08-21 01:12:05 +05:30
|
|
|
if filenames == []: filenames = [None]
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2015-08-08 00:33:26 +05:30
|
|
|
for name in filenames:
|
2015-10-13 22:32:07 +05:30
|
|
|
damask.util.report(scriptName,name)
|
2015-08-08 00:33:26 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
|
2015-08-08 00:33:26 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
randomSeed = int(os.urandom(4).hex(),16) if options.randomSeed is None else options.randomSeed # random seed per file
|
|
|
|
random.seed(randomSeed)
|
2015-04-27 11:00:29 +05:30
|
|
|
|
|
|
|
# --------------- figure out limits (left/right), delta, and interval -----------------------------
|
|
|
|
ODF = {}
|
2020-03-18 18:17:09 +05:30
|
|
|
eulers = table.get('euler')
|
|
|
|
limits = np.array([np.min(eulers,axis=0),np.max(eulers,axis=0)]) # min/max euler angles in degrees
|
2015-11-12 01:21:34 +05:30
|
|
|
ODF['limit'] = np.radians(limits[1,:]) # right hand limits in radians
|
2015-08-21 01:12:05 +05:30
|
|
|
ODF['center'] = 0.0 if all(limits[0,:]<1e-8) else 0.5 # vertex or cell centered
|
2015-04-27 11:00:29 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
ODF['interval'] = np.array(list(map(len,[np.unique(eulers[:,i]) for i in range(3)])),'i') # steps are number of distict values
|
2015-04-27 11:00:29 +05:30
|
|
|
ODF['nBins'] = ODF['interval'].prod()
|
2015-11-12 01:21:34 +05:30
|
|
|
ODF['delta'] = np.radians(np.array(limits[1,0:3]-limits[0,0:3])/(ODF['interval']-1)) # step size
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
if eulers.shape[0] != ODF['nBins']:
|
|
|
|
damask.util.croak('expecting %i values but got %i'%(ODF['nBins'],eulers.shape[0]))
|
2015-04-27 11:00:29 +05:30
|
|
|
continue
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2015-11-12 01:21:34 +05:30
|
|
|
# ----- build binnedODF array and normalize ------------------------------------------------------
|
2015-04-15 20:07:46 +05:30
|
|
|
sumdV_V = 0.0
|
|
|
|
ODF['dV_V'] = [None]*ODF['nBins']
|
|
|
|
ODF['nNonZero'] = 0
|
2015-04-27 11:00:29 +05:30
|
|
|
dg = ODF['delta'][0]*2.0*math.sin(ODF['delta'][1]/2.0)*ODF['delta'][2]
|
2020-03-18 18:17:09 +05:30
|
|
|
intensity = table.get('intensity')
|
2015-04-27 11:00:29 +05:30
|
|
|
for b in range(ODF['nBins']):
|
|
|
|
ODF['dV_V'][b] = \
|
2020-03-18 18:17:09 +05:30
|
|
|
max(0.0,intensity[b,0]) * dg * \
|
2015-08-08 00:33:26 +05:30
|
|
|
math.sin(((b//ODF['interval'][2])%ODF['interval'][1]+ODF['center'])*ODF['delta'][1])
|
2015-04-27 11:00:29 +05:30
|
|
|
if ODF['dV_V'][b] > 0.0:
|
|
|
|
sumdV_V += ODF['dV_V'][b]
|
2015-04-15 20:07:46 +05:30
|
|
|
ODF['nNonZero'] += 1
|
|
|
|
|
2015-10-13 22:32:07 +05:30
|
|
|
for b in range(ODF['nBins']):
|
|
|
|
ODF['dV_V'][b] /= sumdV_V # normalize dV/V
|
2015-08-08 00:33:26 +05:30
|
|
|
|
2015-10-13 22:32:07 +05:30
|
|
|
damask.util.croak(['non-zero fraction: %12.11f (%i/%i)'%(float(ODF['nNonZero'])/ODF['nBins'],
|
2015-08-08 00:33:26 +05:30
|
|
|
ODF['nNonZero'],
|
|
|
|
ODF['nBins']),
|
|
|
|
'Volume integral of ODF: %12.11f\n'%sumdV_V,
|
|
|
|
'Reference Integral: %12.11f\n'%(ODF['limit'][0]*ODF['limit'][2]*(1-math.cos(ODF['limit'][1]))),
|
|
|
|
])
|
|
|
|
|
2015-04-15 20:07:46 +05:30
|
|
|
Functions = {'IA': 'directInversion', 'STAT': 'TothVanHoutteSTAT', 'MC': 'MonteCarloBins'}
|
|
|
|
method = Functions[options.algorithm]
|
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
Orientations, ReconstructedODF = (globals()[method])(ODF,options.number)
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2016-03-02 17:13:09 +05:30
|
|
|
# calculate accuracy of sample
|
2015-04-15 20:07:46 +05:30
|
|
|
squaredDiff = {'orig':0.0,method:0.0}
|
|
|
|
squaredRelDiff = {'orig':0.0,method:0.0}
|
|
|
|
mutualProd = {'orig':0.0,method:0.0}
|
|
|
|
indivSum = {'orig':0.0,method:0.0}
|
|
|
|
indivSquaredSum = {'orig':0.0,method:0.0}
|
|
|
|
|
|
|
|
for bin in range(ODF['nBins']):
|
|
|
|
squaredDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[bin])**2
|
|
|
|
if ODF['dV_V'][bin] > 0.0:
|
|
|
|
squaredRelDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[bin])**2/ODF['dV_V'][bin]**2
|
|
|
|
mutualProd[method] += ODF['dV_V'][bin]*ReconstructedODF[bin]
|
|
|
|
indivSum[method] += ReconstructedODF[bin]
|
|
|
|
indivSquaredSum[method] += ReconstructedODF[bin]**2
|
2015-03-05 15:35:00 +05:30
|
|
|
indivSum['orig'] += ODF['dV_V'][bin]
|
|
|
|
indivSquaredSum['orig'] += ODF['dV_V'][bin]**2
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
damask.util.croak(['sqrt(N*)RMSD of ODFs:\t %12.11f'% math.sqrt(options.number*squaredDiff[method]),
|
2015-08-08 00:33:26 +05:30
|
|
|
'RMSrD of ODFs:\t %12.11f'%math.sqrt(squaredRelDiff[method]),
|
|
|
|
'rMSD of ODFs:\t %12.11f'%(squaredDiff[method]/indivSquaredSum['orig']),
|
|
|
|
'nNonZero correlation slope:\t %12.11f'\
|
2015-04-27 11:00:29 +05:30
|
|
|
%((ODF['nNonZero']*mutualProd[method]-indivSum['orig']*indivSum[method])/\
|
2015-08-08 00:33:26 +05:30
|
|
|
(ODF['nNonZero']*indivSquaredSum['orig']-indivSum['orig']**2)),
|
|
|
|
'nNonZero correlation confidence:\t %12.11f'\
|
2015-04-27 11:00:29 +05:30
|
|
|
%((mutualProd[method]-indivSum['orig']*indivSum[method]/ODF['nNonZero'])/\
|
2015-08-08 00:33:26 +05:30
|
|
|
(ODF['nNonZero']*math.sqrt((indivSquaredSum['orig']/ODF['nNonZero']-(indivSum['orig']/ODF['nNonZero'])**2)*\
|
|
|
|
(indivSquaredSum[method]/ODF['nNonZero']-(indivSum[method]/ODF['nNonZero'])**2)))),
|
|
|
|
])
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
if method == 'IA' and options.number < ODF['nNonZero']:
|
2015-03-05 15:35:00 +05:30
|
|
|
strOpt = '(%i)'%ODF['nNonZero']
|
2015-03-06 00:31:31 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
formatwidth = 1+int(math.log10(options.number))
|
2015-08-08 00:33:26 +05:30
|
|
|
|
|
|
|
materialConfig = [
|
|
|
|
'#' + scriptID + ' ' + ' '.join(sys.argv[1:]),
|
2015-10-13 22:32:07 +05:30
|
|
|
'# random seed %i'%randomSeed,
|
2015-08-08 00:33:26 +05:30
|
|
|
'#-------------------#',
|
|
|
|
'<microstructure>',
|
|
|
|
'#-------------------#',
|
|
|
|
]
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
for i,ID in enumerate(range(options.number)):
|
2015-08-08 00:33:26 +05:30
|
|
|
materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)),
|
|
|
|
'(constituent) phase %i texture %s fraction 1.0'%(options.phase,str(ID+1).rjust(formatwidth)),
|
|
|
|
]
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2015-08-08 00:33:26 +05:30
|
|
|
materialConfig += [
|
|
|
|
'#-------------------#',
|
|
|
|
'<texture>',
|
|
|
|
'#-------------------#',
|
|
|
|
]
|
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
for ID in range(options.number):
|
2015-08-08 00:33:26 +05:30
|
|
|
eulers = Orientations[ID]
|
2015-04-15 20:07:46 +05:30
|
|
|
|
2015-08-08 00:33:26 +05:30
|
|
|
materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)),
|
2015-10-13 22:32:07 +05:30
|
|
|
'(gauss) phi1 {} Phi {} phi2 {} scatter 0.0 fraction 1.0'.format(*eulers),
|
2015-08-08 00:33:26 +05:30
|
|
|
]
|
|
|
|
|
2016-03-02 17:13:09 +05:30
|
|
|
#--- output finalization --------------------------------------------------------------------------
|
2015-08-08 00:33:26 +05:30
|
|
|
|
2020-03-18 18:17:09 +05:30
|
|
|
with (open(os.path.splitext(name)[0]+'_'+method+'_'+str(options.number)+'_material.config','w')) as outfile:
|
2015-08-08 00:33:26 +05:30
|
|
|
outfile.write('\n'.join(materialConfig)+'\n')
|