Merge branch 'orientationClass_with_negative_P' into development

This commit is contained in:
Martin Diehl 2018-12-20 20:43:00 +01:00
commit 1ef82e351b
17 changed files with 349 additions and 609 deletions

@ -1 +1 @@
Subproject commit 58137906b84b6cf0e273dfdde623a2986d03f98e
Subproject commit b9a52a85cd65cc27a8e863302bd984abdcad1455

View File

@ -13,7 +13,7 @@ from .asciitable import ASCIItable # noqa
from .config import Material # noqa
from .colormaps import Colormap, Color # noqa
from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
from .orientation import Quaternion, Symmetry, Orientation # noqa
#from .block import Block # only one class
from .result import Result # noqa

View File

@ -494,7 +494,7 @@ class ASCIItable():
1)))
use = np.array(columns) if len(columns) > 0 else None
self.tags = list(np.array(self.tags)[use]) # update labels with valid subset
self.tags = list(np.array(self.__IO__['tags'])[use]) # update labels with valid subset
self.data = np.loadtxt(self.__IO__['in'],usecols=use,ndmin=2)
# self.data = np.genfromtxt(self.__IO__['in'],dtype=None,names=self.tags,usecols=use)

View File

@ -7,24 +7,6 @@
import math,os
import numpy as np
# ******************************************************************************************
class Rodrigues:
def __init__(self, vector = np.zeros(3)):
self.vector = vector
def asQuaternion(self):
norm = np.linalg.norm(self.vector)
halfAngle = np.arctan(norm)
return Quaternion(np.cos(halfAngle),np.sin(halfAngle)*self.vector/norm)
def asAngleAxis(self):
norm = np.linalg.norm(self.vector)
halfAngle = np.arctan(norm)
return (2.0*halfAngle,self.vector/norm)
# ******************************************************************************************
class Quaternion:
u"""
@ -87,11 +69,11 @@ class Quaternion:
"""Multiplication"""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
try: # quaternion
try: # quaternion
return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
except: pass
try: # vector (perform passive rotation)
try: # vector (perform passive rotation)
( x, y, z) = self.p
(Vx,Vy,Vz) = other[0:3]
A = self.q*self.q - np.dot(self.p,self.p)
@ -104,7 +86,7 @@ class Quaternion:
A*Vz + B*z + C*(x*Vy - y*Vx),
])
except: pass
try: # scalar
try: # scalar
return self.__class__(q=self.q*other,
p=self.p*other)
except:
@ -114,7 +96,7 @@ class Quaternion:
"""In-place multiplication"""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
try: # Quaternion
try: # Quaternion
self.q = self.q*other.q - np.dot(self.p,other.p)
self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
except: pass
@ -178,12 +160,13 @@ class Quaternion:
magnitude = __abs__
def __eq__(self,other):
"""Equal at e-8 precision"""
return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8
"""Equal (sufficiently close) to each other"""
return np.isclose(( self-other).magnitude(),0.0) \
or np.isclose((-self-other).magnitude(),0.0)
def __ne__(self,other):
"""Not equal at e-8 precision"""
return not self.__eq__(self,other)
"""Not equal (sufficiently close) to each other"""
return not self.__eq__(other)
def __cmp__(self,other):
"""Linear ordering"""
@ -193,11 +176,6 @@ class Quaternion:
def magnitude_squared(self):
return self.q ** 2 + np.dot(self.p,self.p)
def identity(self):
self.q = 1.
self.p = np.zeros(3,dtype=float)
return self
def normalize(self):
d = self.magnitude()
if d > 0.0:
@ -209,13 +187,6 @@ class Quaternion:
self.p = -self.p
return self
def inverse(self):
d = self.magnitude()
if d > 0.0:
self.conjugate()
self /= d
return self
def homomorph(self):
if self.q < 0.0:
self.q = -self.q
@ -228,16 +199,13 @@ class Quaternion:
def conjugated(self):
return self.copy().conjugate()
def inversed(self):
return self.copy().inverse()
def homomorphed(self):
return self.copy().homomorph()
def asList(self):
return [self.q]+list(self.p)
def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.)
def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.)
return np.outer(self.asList(),self.asList())
def asMatrix(self):
@ -257,24 +225,26 @@ class Quaternion:
])
def asAngleAxis(self,
degrees = False):
if self.q > 1.:
self.normalize()
degrees = False,
flat = False):
s = math.sqrt(1. - self.q**2)
x = 2*self.q**2 - 1.
y = 2*self.q * s
angle = 2.0*math.acos(self.q)
angle = math.atan2(y,x)
if angle < 0.0:
angle *= -1.
s *= -1.
if np.isclose(angle,0.0):
angle = 0.0
axis = np.array([0.0,0.0,1.0])
elif np.isclose(self.q,0.0):
angle = math.pi
axis = self.p
else:
axis = np.sign(self.q)*self.p/np.linalg.norm(self.p)
return (np.degrees(angle) if degrees else angle,
np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))
angle = np.degrees(angle) if degrees else angle
return np.hstack((angle,axis)) if flat else (angle,axis)
def asRodrigues(self):
return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q
return np.inf*np.ones(3) if np.isclose(self.q,0.0) else self.p/self.q
def asEulers(self,
degrees = False):
@ -285,9 +255,9 @@ class Quaternion:
q12 = self.p[0]**2 + self.p[1]**2
chi = np.sqrt(q03*q12)
if abs(chi) < 1e-10 and abs(q12) < 1e-10:
if np.isclose(chi,0.0) and np.isclose(q12,0.0):
eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0])
elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
elif np.isclose(chi,0.0) and np.isclose(q03,0.0):
eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0])
else:
eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi),
@ -295,6 +265,7 @@ class Quaternion:
math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi),
])
eulers %= 2.0*math.pi # enforce positive angles
return np.degrees(eulers) if degrees else eulers
@ -311,10 +282,12 @@ class Quaternion:
randomSeed = int(binascii.hexlify(os.urandom(4)),16)
np.random.seed(randomSeed)
r = np.random.random(3)
w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2])
x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2])
A = math.sqrt(max(0.0,r[2]))
B = math.sqrt(max(0.0,1.0-r[2]))
w = math.cos(2.0*math.pi*r[0])*A
x = math.sin(2.0*math.pi*r[1])*B
y = math.cos(2.0*math.pi*r[1])*B
z = math.sin(2.0*math.pi*r[0])*A
return cls(quat=[w,x,y,z])
@ -372,10 +345,10 @@ class Quaternion:
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
y = P*0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2])
z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2])
w = 0.5*math.sqrt(max(0.0,1.0+m[0,0]+m[1,1]+m[2,2]))
x = P*0.5*math.sqrt(max(0.0,1.0+m[0,0]-m[1,1]-m[2,2]))
y = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]+m[1,1]-m[2,2]))
z = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]-m[1,1]+m[2,2]))
x *= -1 if m[2,1] < m[1,2] else 1
y *= -1 if m[0,2] < m[2,0] else 1
@ -443,16 +416,16 @@ class Symmetry:
def __repr__(self):
"""Readbable string"""
"""Readable string"""
return '{}'.format(self.lattice)
def __eq__(self, other):
"""Equal"""
"""Equal to other"""
return self.lattice == other.lattice
def __neq__(self, other):
"""Not equal"""
"""Not equal to other"""
return not self.__eq__(other)
def __cmp__(self,other):
@ -529,7 +502,7 @@ class Symmetry:
]
return list(map(Quaternion,
np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))]))
np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))]))
def equivalentQuaternions(self,
@ -541,7 +514,7 @@ class Symmetry:
def inFZ(self,R):
"""Check whether given Rodrigues vector falls into fundamental zone of own symmetry."""
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentally passed quaternion
# fundamental zone in Rodrigues space is point symmetric around origin
R = abs(R)
if self.lattice == 'cubic':
@ -668,7 +641,7 @@ class Symmetry:
if color: # have to return color array
if inSST:
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps
rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity
rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1
else:
rgb = np.zeros(3,dtype=float)
@ -747,8 +720,9 @@ class Orientation:
rodrigues = property(asRodrigues)
def asAngleAxis(self,
degrees = False):
return self.quaternion.asAngleAxis(degrees)
degrees = False,
flat = False):
return self.quaternion.asAngleAxis(degrees,flat)
angleAxis = property(asAngleAxis)
def asMatrix(self):

View File

@ -132,6 +132,30 @@ class extendableOption(Option):
else:
Option.take_action(self, action, dest, opt, value, values, parser)
# Print iterations progress
# from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
def progressBar(iteration, total, prefix='', suffix='', decimals=1, bar_length=100):
"""
Call in a loop to create terminal progress bar
@params:
iteration - Required : current iteration (Int)
total - Required : total iterations (Int)
prefix - Optional : prefix string (Str)
suffix - Optional : suffix string (Str)
decimals - Optional : positive number of decimals in percent complete (Int)
bar_length - Optional : character length of bar (Int)
"""
str_format = "{0:." + str(decimals) + "f}"
percents = str_format.format(100 * (iteration / float(total)))
filled_length = int(round(bar_length * iteration / float(total)))
bar = '' * filled_length + '-' * (bar_length - filled_length)
sys.stderr.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)),
if iteration == total: sys.stderr.write('\n\n')
sys.stderr.flush()
# -----------------------------
class backgroundMessage(threading.Thread):
"""Reporting with animation to indicate progress"""

View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
@ -19,55 +19,50 @@ Transform X,Y,Z,F APS BeamLine 34 coordinates to x,y,z APS strain coordinates.
""", version = scriptID)
parser.add_option('-f','--frame', dest='frame', nargs=4, type='string', metavar='string string string string',
help='APS X,Y,Z coords, and depth F')
parser.set_defaults(frame = None)
parser.add_option('-f',
'--frame',
dest='frame',
metavar='string',
help='APS X,Y,Z coords')
parser.add_option('--depth',
dest='depth',
metavar='string',
help='depth')
(options,filenames) = parser.parse_args()
if options.frame is None:
parser.error('no data column specified...')
parser.error('frame not specified')
if options.depth is None:
parser.error('depth not specified')
# --- loop over input files ------------------------------------------------------------------------
datainfo = {'len':3,
'label':[]
}
datainfo['label'] += options.frame
# --- loop over input files -------------------------------------------------------------------------
if filenames == []:
filenames = ['STDIN']
if filenames == []: filenames = [None]
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+'_tmp','w'), 'croak':sys.stderr}
file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n')
try: table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks -----------------------------------------
errors = []
if table.label_dimension(options.frame) != 3:
errors.append('input {} does not have dimension 3.'.format(options.frame))
if table.label_dimension(options.depth) != 1:
errors.append('input {} does not have dimension 1.'.format(options.depth))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
table = damask.ASCIItable(file['input'],file['output'],buffered=False) # make unbuffered ASCII_table
table.head_read() # read ASCII header info
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
# --------------- figure out columns to process ---------------------------------------------------
active = []
column = {}
columnMissing = False
for label in datainfo['label']:
key = label
if key in table.labels(raw = True):
active.append(label)
column[label] = table.labels.index(key) # remember columns of requested data
else:
file['croak'].write('column %s not found...\n'%label)
columnMissing = True
if columnMissing: continue
# ------------------------------------------ assemble header ---------------------------------------
table.labels_append(['%i_coord'%(i+1) for i in range(3)]) # extend ASCII header with new labels
table.head_write()
@ -77,21 +72,15 @@ for name in filenames:
RotMat2TSL=np.array([[1., 0., 0.],
[0., np.cos(theta), np.sin(theta)], # Orientation to account for -135 deg
[0., -np.sin(theta), np.cos(theta)]]) # rotation for TSL convention
vec = np.zeros(4)
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for i,label in enumerate(active):
vec[i] = table.data[column[label]]
coord = list(map(float,table.data[table.label_index(options.frame):table.label_index(options.frame)+3]))
depth = float(table.data[table.label_index(options.depth)])
table.data_append(np.dot(RotMat2TSL,np.array([-vec[0], -vec[1],-vec[2]+vec[3]])))
table.data_append(np.dot(RotMat2TSL,np.array([-coord[0],-coord[1],-coord[2]+depth])))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output result -----------------------------------------
outputAlive and table.output_flush() # just in case of buffered ASCII table
# ------------------------------------------ output finalization -----------------------------------
table.input_close() # close input ASCII table (works for stdin)
table.output_close() # close output ASCII table (works for stdout)
if file['name'] != 'STDIN':
os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new
table.close() # close ASCII tables

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,time,copy
import os,sys,copy
import numpy as np
import damask
from optparse import OptionParser
@ -29,49 +29,28 @@ parser.add_option('-d',
parser.add_option('-s',
'--symmetry',
dest = 'symmetry',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'crystal symmetry [%default]')
parser.add_option('-e',
'--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'label of Euler angles')
parser.add_option('--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m',
'--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'label of orientation matrix')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'label of crystal frame a vector')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'label of crystal frame b vector')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'label of crystal frame c vector')
parser.add_option('-q',
'--quaternion',
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'label of quaternion')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.add_option('-p',
'--pos', '--position',
dest = 'pos',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('--quiet',
dest='verbose',
action = 'store_false',
help = 'hide status bar (useful when piping to file)')
parser.set_defaults(disorientation = 5,
verbose = True,
quaternion = 'orientation',
symmetry = 'cubic',
pos = 'pos',
degrees = False,
)
(options, filenames) = parser.parse_args()
@ -79,22 +58,6 @@ parser.set_defaults(disorientation = 5,
if options.radius is None:
parser.error('no radius specified.')
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
cos_disorientation = np.cos(np.radians(options.disorientation/2.)) # cos of half the disorientation angle
# --- loop over input files -------------------------------------------------------------------------
@ -118,9 +81,9 @@ for name in filenames:
if not 3 >= table.label_dimension(options.pos) >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
if not np.all(table.label_dimension(label) == dim):
errors.append('input "{}" does not have dimension {}.'.format(label,dim))
else: column = table.label_index(label)
if not np.all(table.label_dimension(options.quaternion) == 4):
errors.append('input "{}" does not have dimension 4.'.format(options.quaternion))
else: column = table.label_index(options.quaternion)
if remarks != []: damask.util.croak(remarks)
if errors != []:
@ -131,34 +94,18 @@ for name in filenames:
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append('grainID_{}@{:g}'.format('+'.join(label)
if isinstance(label, (list,tuple))
else label,
options.disorientation)) # report orientation source and disorientation
table.labels_append('grainID_{}@{:g}'.format(options.quaternion,options.disorientation)) # report orientation source and disorientation
table.head_write()
# ------------------------------------------ process data ------------------------------------------
# ------------------------------------------ build KD tree -----------------------------------------
# --- start background messaging
bg = damask.util.backgroundMessage()
bg.start()
bg.set_message('reading positions...')
table.data_readArray(options.pos) # read position vectors
grainID = -np.ones(len(table.data),dtype=int)
start = tick = time.clock()
bg.set_message('building KD tree...')
Npoints = table.data.shape[0]
kdtree = spatial.KDTree(copy.deepcopy(table.data))
# ------------------------------------------ assign grain IDs --------------------------------------
tick = time.clock()
orientations = [] # quaternions found for grain
memberCounts = [] # number of voxels in grain
p = 0 # point counter
@ -169,26 +116,11 @@ for name in filenames:
table.data_rewind()
while table.data_read(): # read next data line of ASCII table
if p > 0 and p % 1000 == 0:
if options.verbose and Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=p,total=Npoints)
time_delta = (time.clock()-tick) * (len(grainID) - p) / p
bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'\
%(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),np.count_nonzero(memberCounts)))
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians,
symmetry = options.symmetry).reduced()
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3])).reshape(3,3),
symmetry = options.symmetry).reduced()
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])),
symmetry = options.symmetry).reduced()
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
matched = False
alreadyChecked = {}
@ -233,13 +165,12 @@ for name in filenames:
outputAlive = True
p = 0
damask.util.progressBar(iteration=1,total=1)
while outputAlive and table.data_read(): # read next data line of ASCII table
table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID
outputAlive = table.data_write() # output processed line
p += 1
bg.set_message('done after {} seconds'.format(time.clock()-start))
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os,sys
import numpy as np
from optparse import OptionParser
import damask
@ -18,66 +18,29 @@ Add RGB color value corresponding to TSL-OIM scheme for inverse pole figures.
""", version = scriptID)
parser.add_option('-p', '--pole',
parser.add_option('-p',
'--pole',
dest = 'pole',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'lab frame direction for inverse pole figure [%default]')
parser.add_option('-s', '--symmetry',
parser.add_option('-s',
'--symmetry',
dest = 'symmetry',
type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string',
help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion',
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'quaternion label')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(pole = (0.0,0.0,1.0),
quaternion = 'orientation',
symmetry = damask.Symmetry.lattices[-1],
degrees = False,
)
(options, filenames) = parser.parse_args()
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
pole = np.array(options.pole)
pole /= np.linalg.norm(pole)
@ -98,12 +61,12 @@ for name in filenames:
# ------------------------------------------ sanity checks ----------------------------------------
if not np.all(table.label_dimension(label) == dim):
damask.util.croak('input {} does not have dimension {}.'.format(label,dim))
if not table.label_dimension(options.quaternion) == 4:
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
table.close(dismiss = True) # close ASCIItable and remove empty file
continue
column = table.label_index(label)
column = table.label_index(options.quaternion)
# ------------------------------------------ assemble header ---------------------------------------
@ -115,20 +78,8 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,
symmetry = options.symmetry).reduced()
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),
symmetry = options.symmetry).reduced()
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
table.data_append(o.IPFcolor(pole))
outputAlive = table.data_write() # output processed line

View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os

View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os,sys
import numpy as np
from optparse import OptionParser
import damask
@ -9,6 +9,31 @@ import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# convention conformity checks
# --------------------------------------------------------------------
def check_Eulers(eulers):
if np.any(eulers < 0.0) or np.any(eulers > 2.0*np.pi) or eulers[1] > np.pi: # Euler angles within valid range?
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers))
return eulers
def check_quaternion(q):
if q[0] < 0.0: # positive first quaternion component?
raise ValueError('quaternion has negative first component.\n{}'.format(q[0]))
if not np.isclose(np.linalg.norm(q), 1.0): # unit quaternion?
raise ValueError('quaternion is not of unit length.\n{} {} {} {}'.format(*q))
return q
def check_matrix(M):
if not np.isclose(np.linalg.det(M),1.0): # proper rotation?
raise ValueError('matrix is not a proper rotation.\n{}'.format(M))
if not np.isclose(np.dot(M[0],M[1]), 0.0) \
or not np.isclose(np.dot(M[1],M[2]), 0.0) \
or not np.isclose(np.dot(M[2],M[0]), 0.0): # all orthogonal?
raise ValueError('matrix is not orthogonal.\n{}'.format(M))
return M
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
@ -21,58 +46,64 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can
""", version = scriptID)
outputChoices = ['quaternion','rodrigues','eulers']
parser.add_option('-o', '--output',
outputChoices = {
'quaternion': ['quat',4],
'rodrigues': ['rodr',3],
'eulers': ['eulr',3],
'matrix': ['mtrx',9],
'angleaxis': ['aaxs',4],
}
parser.add_option('-o',
'--output',
dest = 'output',
action = 'extend', metavar = '<string LIST>',
help = 'output orientation formats {{{}}}'.format(', '.join(outputChoices)))
parser.add_option('-s', '--symmetry',
dest = 'symmetry',
type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string',
help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('-d', '--degrees',
parser.add_option('-d',
'--degrees',
dest = 'degrees',
action = 'store_true',
help = 'angles are given in degrees [%default]')
parser.add_option('-R', '--labrotation',
help = 'all angles in degrees')
parser.add_option('-R',
'--labrotation',
dest='labrotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis of additional lab frame rotation')
parser.add_option('-r', '--crystalrotation',
parser.add_option('-r',
'--crystalrotation',
dest='crystalrotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis of additional crystal frame rotation')
parser.add_option( '--eulers',
parser.add_option('--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'Euler angles label')
parser.add_option( '--rodrigues',
parser.add_option('--rodrigues',
dest = 'rodrigues',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'Rodrigues vector label')
parser.add_option( '--matrix',
parser.add_option('--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'orientation matrix label')
parser.add_option( '--quaternion',
parser.add_option('--quaternion',
dest = 'quaternion',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'quaternion label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-x',
dest = 'x',
metavar = 'string',
help = 'label of lab x vector (expressed in crystal coords)')
parser.add_option('-y',
dest = 'y',
metavar = 'string',
help = 'label of lab y vector (expressed in crystal coords)')
parser.add_option('-z',
dest = 'z',
metavar = 'string',
help = 'label of lab z vector (expressed in crystal coords)')
parser.set_defaults(output = [],
symmetry = damask.Symmetry.lattices[-1],
labrotation = (0.,1.,1.,1.), # no rotation about 1,1,1
crystalrotation = (0.,1.,1.,1.), # no rotation about 1,1,1
degrees = False,
@ -86,9 +117,9 @@ if options.output == [] or (not set(options.output).issubset(set(outputChoices))
input = [options.eulers is not None,
options.rodrigues is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.x is not None and \
options.y is not None and \
options.z is not None,
options.matrix is not None,
options.quaternion is not None,
]
@ -97,13 +128,14 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
(options.rodrigues,3,'rodrigues'),
([options.a,options.b,options.c],[3,3,3],'frame'),
([options.x,options.y,options.z],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
r = damask.Quaternion().fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation
R = damask.Quaternion().fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation
toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation
R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation
# --- loop over input files ------------------------------------------------------------------------
@ -137,32 +169,31 @@ for name in filenames:
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for output in options.output:
if output == 'quaternion': table.labels_append(['{}_{}_{}({})'.format(i+1,'quat',options.symmetry,label) for i in range(4)])
elif output == 'rodrigues': table.labels_append(['{}_{}_{}({})'.format(i+1,'rodr',options.symmetry,label) for i in range(3)])
elif output == 'eulers': table.labels_append(['{}_{}_{}({})'.format(i+1,'eulr',options.symmetry,label) for i in range(3)])
if output in outputChoices:
table.labels_append(['{}_{}({})'.format(i+1,outputChoices[output][0],label) \
for i in range(outputChoices[output][1])])
table.head_write()
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,
symmetry = options.symmetry).reduced()
if inputtype == 'eulers':
o = damask.Orientation(Eulers = check_Eulers(np.array(list(map(float,table.data[column:column+3])))*toRadians))
elif inputtype == 'rodrigues':
o = damask.Orientation(Rodrigues= np.array(list(map(float,table.data[column:column+3]))),
symmetry = options.symmetry).reduced()
o = damask.Orientation(Rodrigues = np.array(list(map(float,table.data[column:column+3]))))
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
o = damask.Orientation(matrix = check_matrix(np.array(list(map(float,table.data[column:column+9]))).reshape(3,3)))
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),
symmetry = options.symmetry).reduced()
M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3).T
o = damask.Orientation(matrix = check_matrix(M/np.linalg.norm(M,axis=0)))
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
o = damask.Orientation(quaternion = check_quaternion(np.array(list(map(float,table.data[column:column+4])))))
o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations
@ -170,6 +201,8 @@ for name in filenames:
if output == 'quaternion': table.data_append(o.asQuaternion())
elif output == 'rodrigues': table.data_append(o.asRodrigues())
elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees))
elif output == 'matrix': table.data_append(o.asMatrix())
elif output == 'angleaxis': table.data_append(o.asAngleAxis(degrees=options.degrees,flat=True))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os,sys
import numpy as np
from optparse import OptionParser
import damask
@ -14,70 +14,32 @@ scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Add x,y coordinates of stereographic projection of given direction (pole) in crystal frame.
Add coordinates of stereographic projection of given direction (pole) in crystal frame.
""", version = scriptID)
parser.add_option('-p', '--pole',
parser.add_option('-p',
'--pole',
dest = 'pole',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'crystal frame direction for pole figure [%default]')
parser.add_option('--polar',
dest = 'polar',
action = 'store_true',
help = 'output polar coordinates r,phi [%default]')
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion',
help = 'output polar coordinates (r,φ) instead of Cartesian coordinates (x,y)')
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'quaternion label')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(pole = (1.0,0.0,0.0),
degrees = False,
quaternion = 'orientation',
polar = False,
)
(options, filenames) = parser.parse_args()
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
pole = np.array(options.pole)
pole /= np.linalg.norm(pole)
@ -98,18 +60,13 @@ for name in filenames:
# ------------------------------------------ sanity checks ----------------------------------------
errors = []
remarks = []
if not np.all(table.label_dimension(label) == dim): errors.append('input {} does not have dimension {}.'.format(label,dim))
else: column = table.label_index(label)
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
if not table.label_dimension(options.quaternion) == 4:
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
table.close(dismiss = True) # close ASCIItable and remove empty file
continue
column = table.label_index(options.quaternion)
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
@ -119,16 +76,7 @@ for name in filenames:
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians)
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose())
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3))
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
rotatedPole = o.quaternion*pole # rotate pole according to crystal orientation
(x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection

View File

@ -109,64 +109,42 @@ Add columns listing Schmid factors (and optional trace vector of selected system
""", version = scriptID)
latticeChoices = ('fcc','bcc','hex')
parser.add_option('-l','--lattice',
parser.add_option('-l',
'--lattice',
dest = 'lattice', type = 'choice', choices = latticeChoices, metavar='string',
help = 'type of lattice structure [%default] {}'.format(latticeChoices))
parser.add_option('--covera',
dest = 'CoverA', type = 'float', metavar = 'float',
help = 'C over A ratio for hexagonal systems')
parser.add_option('-f', '--force',
parser.add_option('-f',
'--force',
dest = 'force',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'force direction in lab frame [%default]')
parser.add_option('-n', '--normal',
parser.add_option('-n',
'--normal',
dest = 'normal',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'stress plane normal in lab frame [%default]')
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion',
help = 'stress plane normal in lab frame, per default perpendicular to the force')
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'quaternion label')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(force = (0.0,0.0,1.0),
quaternion='orientation',
normal = None,
lattice = latticeChoices[0],
CoverA = math.sqrt(8./3.),
degrees = False,
)
(options, filenames) = parser.parse_args()
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
force = np.array(options.force)
force /= np.linalg.norm(force)
if options.normal:
if options.normal is not None:
normal = np.array(options.normal)
normal /= np.linalg.norm(normal)
if abs(np.dot(force,normal)) > 1e-3:
@ -174,22 +152,6 @@ if options.normal:
else:
normal = force
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
slip_normal = np.zeros_like(slip_direction)
@ -227,13 +189,12 @@ for name in filenames:
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
if not np.all(table.label_dimension(label) == dim):
damask.util.croak('input {} does not have dimension {}.'.format(label,dim))
if not table.label_dimension(options.quaternion) == 4:
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
table.close(dismiss = True) # close ASCIItable and remove empty file
continue
column = table.label_index(label)
column = table.label_index(options.quaternion)
# ------------------------------------------ assemble header ---------------------------------------
@ -251,17 +212,7 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,)
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),)
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),)
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),)
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
table.data_append( np.abs( np.sum(slip_direction * (o.quaternion * force) ,axis=1) \
* np.sum(slip_normal * (o.quaternion * normal),axis=1)))

View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
@ -118,10 +118,9 @@ for name in filenames:
minmax[c] = np.log(minmax[c]) # change minmax to log, too
delta = minmax[:,1]-minmax[:,0]
(grid,xedges,yedges) = np.histogram2d(table.data[:,0],table.data[:,1],
bins=options.bins,
range=minmax,
range=minmax[0:2,0:2],
weights=None if options.weight is None else table.data[:,2])
if options.normCol:

View File

@ -121,12 +121,8 @@ class MPIEspectral_result: # mimic py_post result object
self._logscales = self._keyedPackedArray('logscales',count=self.N_loadcases,type='i')
self.size = self._keyedPackedArray('size:',count=3,type='d')
if self.size == [None,None,None]: # no 'size' found, try legacy alias 'dimension'
self.size = self._keyedPackedArray('dimension',count=3,type='d')
self.grid = self._keyedPackedArray('grid:',count=3,type='i')
if self.grid == [None,None,None]: # no 'grid' found, try legacy alias 'resolution'
self.grid = self._keyedPackedArray('resolution',count=3,type='i')
self.N_nodes = (self.grid[0]+1)*(self.grid[1]+1)*(self.grid[2]+1)
self.N_elements = self.grid[0] * self.grid[1] * self.grid[2]
@ -142,13 +138,8 @@ class MPIEspectral_result: # mimic py_post result object
# parameters for file handling depending on output format
if options.legacy:
self.tagLen=8
self.fourByteLimit = 2**31 -1 -8
else:
self.tagLen=0
self.tagLen=0
self.expectedFileSize = self.dataOffset+self.N_increments*(self.tagLen+self.N_elements*self.N_element_scalars*8)
if options.legacy: self.expectedFileSize+=self.expectedFileSize//self.fourByteLimit*8 # add extra 8 bytes for additional headers at 4 GB limits
if self.expectedFileSize != self.filesize:
print('\n**\n* Unexpected file size. Incomplete simulation or file corrupted!\n**')
@ -280,42 +271,16 @@ class MPIEspectral_result: # mimic py_post result object
return self.N_element_scalars
def element_scalar(self,e,idx):
if not options.legacy:
incStart = self.dataOffset \
+ self.position*8*self.N_elements*self.N_element_scalars
where = (e*self.N_element_scalars + idx)*8
try:
self.file.seek(incStart+where)
value = struct.unpack('d',self.file.read(8))[0]
except:
print('seeking {}'.format(incStart+where))
print('e {} idx {}'.format(e,idx))
sys.exit(1)
else:
self.fourByteLimit = 2**31 -1 -8
# header & footer + extra header and footer for 4 byte int range (Fortran)
# values
incStart = self.dataOffset \
+ self.position*8*( 1 + self.N_elements*self.N_element_scalars*8//self.fourByteLimit \
+ self.N_elements*self.N_element_scalars)
where = (e*self.N_element_scalars + idx)*8
try:
if where%self.fourByteLimit + 8 >= self.fourByteLimit: # danger of reading into fortran record footer at 4 byte limit
data=''
for i in range(8):
self.file.seek(incStart+where+(where//self.fourByteLimit)*8+4)
data += self.file.read(1)
where += 1
value = struct.unpack('d',data)[0]
else:
self.file.seek(incStart+where+(where//self.fourByteLimit)*8+4)
value = struct.unpack('d',self.file.read(8))[0]
except:
print('seeking {}'.format(incStart+where+(where//self.fourByteLimit)*8+4))
print('e {} idx {}'.format(e,idx))
sys.exit(1)
incStart = self.dataOffset \
+ self.position*8*self.N_elements*self.N_element_scalars
where = (e*self.N_element_scalars + idx)*8
try:
self.file.seek(incStart+where)
value = struct.unpack('d',self.file.read(8))[0]
except:
print('seeking {}'.format(incStart+where))
print('e {} idx {}'.format(e,idx))
sys.exit(1)
return [elemental_scalar(node,value) for node in self.element(e).items]
@ -645,8 +610,6 @@ of already processed data points for evaluation.
parser.add_option('-i','--info', action='store_true', dest='info',
help='list contents of resultfile')
parser.add_option('-l','--legacy', action='store_true', dest='legacy',
help='data format of spectral solver is in legacy format (no MPI out)')
parser.add_option('-n','--nodal', action='store_true', dest='nodal',
help='data is extrapolated to nodal value')
parser.add_option( '--prefix', dest='prefix',
@ -673,10 +636,7 @@ parser.add_option('-p','--type', dest='filetype',
help = 'type of result file [auto]')
parser.add_option('-q','--quiet', dest='verbose',
action = 'store_false',
help = 'suppress verbose output')
parser.add_option('--verbose', dest='verbose',
action = 'store_true',
help = 'enable verbose output')
help = 'hide status bar (useful when piping to file)')
group_material = OptionGroup(parser,'Material identifier')
@ -718,9 +678,8 @@ parser.add_option_group(group_general)
parser.add_option_group(group_special)
parser.set_defaults(info = False,
verbose = False,
legacy = False,
nodal = False,
verbose = True,
prefix = '',
suffix = '',
dir = 'postProc',
@ -747,6 +706,8 @@ if files == []:
parser.print_help()
parser.error('no file specified...')
damask.util.report(scriptName,files[0])
if not os.path.exists(files[0]):
parser.print_help()
parser.error('invalid file "%s" specified...'%files[0])
@ -803,12 +764,6 @@ if not options.constitutiveResult: options.constitutiveResult = []
options.sort.reverse()
options.sep.reverse()
# --- start background messaging
if options.verbose:
bg = damask.util.backgroundMessage()
bg.start()
# --- parse .output and .t16 files
if os.path.splitext(files[0])[1] == '':
@ -825,18 +780,13 @@ me = {
'Constitutive': options.phase,
}
if options.verbose: bg.set_message('parsing .output files...')
for what in me:
outputFormat[what] = ParseOutputFormat(filename, what, me[what])
if '_id' not in outputFormat[what]['specials']:
print("\nsection '{}' not found in <{}>".format(me[what], what))
print('\n'.join(map(lambda x:' [%s]'%x, outputFormat[what]['specials']['brothers'])))
if options.verbose: bg.set_message('opening result file...')
p = OpenPostfile(filename+extension,options.filetype,options.nodal)
if options.verbose: bg.set_message('parsing result file...')
stat = ParsePostfile(p, filename, outputFormat)
if options.filetype == 'marc':
stat['NumberOfIncrements'] -= 1 # t16 contains one "virtual" increment (at 0)
@ -879,8 +829,10 @@ if options.info:
# --- build connectivity maps
elementsOfNode = {}
for e in range(stat['NumberOfElements']):
if options.verbose and e%1000 == 0: bg.set_message('connect elem %i...'%e)
Nelems = stat['NumberOfElements']
for e in range(Nelems):
if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=e,total=Nelems,prefix='1/3: connecting elements')
for n in map(p.node_sequence,p.element(e).items):
if n not in elementsOfNode:
elementsOfNode[n] = [p.element_id(e)]
@ -899,10 +851,13 @@ index = {}
groups = []
groupCount = 0
memberCount = 0
damask.util.progressBar(iteration=1,total=1,prefix='1/3: connecting elements')
if options.nodalScalar:
for n in range(stat['NumberOfNodes']):
if options.verbose and n%1000 == 0: bg.set_message('scan node %i...'%n)
Npoints = stat['NumberOfNodes']
for n in range(Npoints):
if options.verbose and Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=n,total=Npoints,prefix='2/3: scanning nodes ')
myNodeID = p.node_id(n)
myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z]
myElemID = 0
@ -911,32 +866,35 @@ if options.nodalScalar:
# generate an expression that is only true for the locations specified by options.filter
filter = substituteLocation(options.filter, [myElemID,myNodeID,myIpID,myGrainID], myNodeCoordinates)
if filter != '' and not eval(filter): # for all filter expressions that are not true:...
continue # ... ignore this data point and continue with next
if filter != '' and not eval(filter): # for all filter expressions that are not true:...
continue # ... ignore this data point and continue with next
# --- group data locations
# generate a unique key for a group of separated data based on the separation criterium for the location
grp = substituteLocation('#'.join(options.sep), [myElemID,myNodeID,myIpID,myGrainID], myNodeCoordinates)
if grp not in index: # create a new group if not yet present
if grp not in index: # create a new group if not yet present
index[grp] = groupCount
groups.append([[0,0,0,0,0.0,0.0,0.0]]) # initialize with avg location
groups.append([[0,0,0,0,0.0,0.0,0.0]]) # initialize with avg location
groupCount += 1
groups[index[grp]][0][:4] = mapIncremental('','unique',
len(groups[index[grp]])-1,
groups[index[grp]][0][:4],
[myElemID,myNodeID,myIpID,myGrainID]) # keep only if unique average location
[myElemID,myNodeID,myIpID,myGrainID]) # keep only if unique average location
groups[index[grp]][0][4:] = mapIncremental('','avg',
len(groups[index[grp]])-1,
groups[index[grp]][0][4:],
myNodeCoordinates) # incrementally update average location
groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member
myNodeCoordinates) # incrementally update average location
groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member
memberCount += 1
damask.util.progressBar(iteration=1,total=1,prefix='2/3: scanning nodes ')
else:
for e in range(stat['NumberOfElements']):
if options.verbose and e%1000 == 0: bg.set_message('scan elem %i...'%e)
Nelems = stat['NumberOfElements']
for e in range(Nelems):
if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=e,total=Nelems,prefix='2/3: scanning elements ')
myElemID = p.element_id(e)
myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z],
list(map(p.node, map(p.node_sequence, p.element(e).items))))))
@ -976,6 +934,7 @@ else:
myIpCoordinates[n]) # incrementally update average location
groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,n]) # append a new list defining each group member
memberCount += 1
damask.util.progressBar(iteration=1,total=1,prefix='2/3: scanning elements ')
# --------------------------- sort groups --------------------------------
@ -1002,7 +961,6 @@ if 'none' not in map(str.lower, options.sort):
theKeys.append('x[0][%i]'%where[criterium])
sortKeys = eval('lambda x:(%s)'%(','.join(theKeys)))
if options.verbose: bg.set_message('sorting groups...')
groups.sort(key = sortKeys) # in-place sorting to save mem
@ -1021,8 +979,6 @@ standard = ['inc'] + \
# --------------------------- loop over positions --------------------------------
if options.verbose: bg.set_message('getting map between positions and increments...')
incAtPosition = {}
positionOfInc = {}
@ -1048,8 +1004,8 @@ increments = [incAtPosition[x] for x in locations] # build list of increments to
time_start = time.time()
Nincs = len([i for i in locations])
for incCount,position in enumerate(locations): # walk through locations
p.moveto(position+offset_pos) # wind to correct position
# --------------------------- file management --------------------------------
@ -1075,16 +1031,14 @@ for incCount,position in enumerate(locations): # walk through locations
# --------------------------- read and map data per group --------------------------------
member = 0
for group in groups:
Ngroups = len(groups)
for j,group in enumerate(groups):
f = incCount*Ngroups + j
if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ')
N = 0 # group member counter
for (e,n,i,g,n_local) in group[1:]: # loop over group members
member += 1
if member%1000 == 0:
time_delta = ((len(locations)*memberCount)/float(member+incCount*memberCount)-1.0)*(time.time()-time_start)
if options.verbose: bg.set_message('(%02i:%02i:%02i) processing point %i of %i from increment %i (position %i)...'
%(time_delta//3600,time_delta%3600//60,time_delta%60,member,memberCount,increments[incCount],position))
newby = [] # current member's data
if options.nodalScalar:
@ -1172,6 +1126,7 @@ for incCount,position in enumerate(locations): # walk through locations
group[0] + \
mappedResult)
)) + '\n')
damask.util.progressBar(iteration=1,total=1,prefix='3/3: processing points ')
if fileOpen:
file.close()

View File

@ -18,19 +18,15 @@ Rotate vector and/or tensor column data by given angle around given axis.
""", version = scriptID)
parser.add_option('-v','--vector',
dest = 'vector',
parser.add_option('-d', '--data',
dest = 'data',
action = 'extend', metavar = '<string LIST>',
help = 'column heading of vector(s) to rotate')
parser.add_option('-t','--tensor',
dest = 'tensor',
action = 'extend', metavar = '<string LIST>',
help = 'column heading of tensor(s) to rotate')
help = 'vector/tensor value(s) label(s)')
parser.add_option('-r', '--rotation',
dest = 'rotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis to rotate data [%default]')
parser.add_option('-d', '--degrees',
parser.add_option('--degrees',
dest = 'degrees',
action = 'store_true',
help = 'angles are given in degrees [%default]')
@ -41,7 +37,7 @@ parser.set_defaults(rotation = (0.,1.,1.,1.),
(options,filenames) = parser.parse_args()
if options.vector is None and options.tensor is None:
if options.data is None:
parser.error('no data column specified.')
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
@ -59,27 +55,24 @@ for name in filenames:
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
items = {
'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []},
'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []},
}
errors = []
remarks = []
column = {}
active = {'vector':[],'tensor':[]}
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
else:
items[type]['active'].append(what)
items[type]['column'].append(table.label_index(what))
for i,dim in enumerate(table.label_dimension(options.data)):
label = options.data[i]
if dim == -1:
remarks.append('"{}" not found...'.format(label))
elif dim == 3:
remarks.append('adding vector "{}"...'.format(label))
active['vector'].append(label)
elif dim == 9:
remarks.append('adding tensor "{}"...'.format(label))
active['tensor'].append(label)
if remarks != []: damask.util.croak(remarks)
if errors != []:
@ -95,20 +88,14 @@ for name in filenames:
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
datatype = 'vector'
for column in items[datatype]['column']: # loop over all requested labels
table.data[column:column+items[datatype]['dim']] = \
q * np.array(list(map(float,table.data[column:column+items[datatype]['dim']])))
datatype = 'tensor'
for column in items[datatype]['column']: # loop over all requested labels
table.data[column:column+items[datatype]['dim']] = \
np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+items[datatype]['dim']]))).\
reshape(items[datatype]['shape']),R.transpose())).reshape(items[datatype]['dim'])
for v in active['vector']:
column = table.label_index(v)
table.data[column:column+3] = q * np.array(list(map(float,table.data[column:column+3])))
for t in active['tensor']:
column = table.label_index(t)
table.data[column:column+9] = \
np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]))).reshape((3,3)),
R.transpose())).reshape((9))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------

View File

@ -60,8 +60,6 @@ eulers = np.array(damask.orientation.Orientation(
degrees = options.degrees,
).asEulers(degrees=True))
damask.util.croak('{} {} {}'.format(*eulers))
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]