Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Martin Diehl 2016-08-08 08:53:55 +02:00
commit 45593b1660
4 changed files with 61 additions and 54 deletions

View File

@ -1 +1 @@
v2.0.1-5-g920cf2c v2.0.1-35-g6c82641

View File

@ -156,11 +156,9 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize field solver information ! initialize field solver information
@ -199,14 +197,14 @@ program DAMASK_spectral
allocate(loadCases(i)%ID(nActiveFields)) allocate(loadCases(i)%ID(nActiveFields))
field = 1 field = 1
loadCases(i)%ID(field) = FIELD_MECH_ID ! mechanical active by default loadCases(i)%ID(field) = FIELD_MECH_ID ! mechanical active by default
if (any(thermal_type == THERMAL_conduction_ID)) then ! thermal field active thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
field = field + 1 field = field + 1
loadCases(i)%ID(field) = FIELD_THERMAL_ID loadCases(i)%ID(field) = FIELD_THERMAL_ID
endif endif thermalActive
if (any(damage_type == DAMAGE_nonlocal_ID)) then ! damage field active damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
field = field + 1 field = field + 1
loadCases(i)%ID(field) = FIELD_DAMAGE_ID loadCases(i)%ID(field) = FIELD_DAMAGE_ID
endif endif damageActive
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

View File

@ -413,11 +413,15 @@ class Quaternion:
@classmethod @classmethod
def fromAngleAxis(cls, angle, axis): def fromAngleAxis(cls,
angle,
axis,
degrees = False):
if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d') if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d')
axis = axis.astype(float)/np.linalg.norm(axis) axis = axis.astype(float)/np.linalg.norm(axis)
s = math.sin(angle / 2.0) angle = np.radians(angle) if degrees else angle
w = math.cos(angle / 2.0) s = math.sin(0.5 * angle)
w = math.cos(0.5 * angle)
x = axis[0] * s x = axis[0] * s
y = axis[1] * s y = axis[1] * s
z = axis[2] * s z = axis[2] * s
@ -425,27 +429,26 @@ class Quaternion:
@classmethod @classmethod
def fromEulers(cls, eulers, type = 'Bunge'): def fromEulers(cls,
eulers,
type = 'Bunge',
degrees = False):
if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d')
eulers = np.radians(eulers) if degrees else eulers
eulers *= 0.5 # reduce to half angles c = np.cos(0.5 * eulers)
s = np.sin(0.5 * eulers)
c1 = math.cos(eulers[0])
s1 = math.sin(eulers[0])
c2 = math.cos(eulers[1])
s2 = math.sin(eulers[1])
c3 = math.cos(eulers[2])
s3 = math.sin(eulers[2])
if type.lower() == 'bunge' or type.lower() == 'zxz': if type.lower() == 'bunge' or type.lower() == 'zxz':
w = c1 * c2 * c3 - s1 * c2 * s3 w = c[0] * c[1] * c[2] - s[0] * c[1] * s[2]
x = c1 * s2 * c3 + s1 * s2 * s3 x = c[0] * s[1] * c[2] + s[0] * s[1] * s[2]
y = - c1 * s2 * s3 + s1 * s2 * c3 y = - c[0] * s[1] * s[2] + s[0] * s[1] * c[2]
z = c1 * c2 * s3 + s1 * c2 * c3 z = c[0] * c[1] * s[2] + s[0] * c[1] * c[2]
else: else:
w = c1 * c2 * c3 - s1 * s2 * s3 w = c[0] * c[1] * c[2] - s[0] * s[1] * s[2]
x = s1 * s2 * c3 + c1 * c2 * s3 x = s[0] * s[1] * c[2] + c[0] * c[1] * s[2]
y = s1 * c2 * c3 + c1 * s2 * s3 y = s[0] * c[1] * c[2] + c[0] * s[1] * s[2]
z = c1 * s2 * c3 - s1 * c2 * s3 z = c[0] * s[1] * c[2] - s[0] * c[1] * s[2]
return cls([w,x,y,z]) return cls([w,x,y,z])
@ -819,6 +822,7 @@ class Orientation:
Eulers = None, Eulers = None,
random = False, # integer to have a fixed seed or True for real random random = False, # integer to have a fixed seed or True for real random
symmetry = None, symmetry = None,
degrees = False,
): ):
if random: # produce random orientation if random: # produce random orientation
if isinstance(random, bool ): if isinstance(random, bool ):
@ -826,16 +830,16 @@ class Orientation:
else: else:
self.quaternion = Quaternion.fromRandom(randomSeed=random) self.quaternion = Quaternion.fromRandom(randomSeed=random)
elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,): # based on given Euler angles elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,): # based on given Euler angles
self.quaternion = Quaternion.fromEulers(Eulers,'bunge') self.quaternion = Quaternion.fromEulers(Eulers,type='bunge',degrees=degrees)
elif isinstance(matrix, np.ndarray) : # based on given rotation matrix elif isinstance(matrix, np.ndarray) : # based on given rotation matrix
self.quaternion = Quaternion.fromMatrix(matrix) self.quaternion = Quaternion.fromMatrix(matrix)
elif isinstance(angleAxis, np.ndarray) and angleAxis.shape == (4,): # based on given angle and rotation axis elif isinstance(angleAxis, np.ndarray) and angleAxis.shape == (4,): # based on given angle and rotation axis
self.quaternion = Quaternion.fromAngleAxis(angleAxis[0],angleAxis[1:4]) self.quaternion = Quaternion.fromAngleAxis(angleAxis[0],angleAxis[1:4],degrees=degrees)
elif isinstance(Rodrigues, np.ndarray) and Rodrigues.shape == (3,): # based on given Rodrigues vector elif isinstance(Rodrigues, np.ndarray) and Rodrigues.shape == (3,): # based on given Rodrigues vector
self.quaternion = Quaternion.fromRodrigues(Rodrigues) self.quaternion = Quaternion.fromRodrigues(Rodrigues)
elif isinstance(quaternion, Quaternion): # based on given quaternion elif isinstance(quaternion, Quaternion): # based on given quaternion
self.quaternion = quaternion.homomorphed() self.quaternion = quaternion.homomorphed()
elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion-like array
self.quaternion = Quaternion(quaternion).homomorphed() self.quaternion = Quaternion(quaternion).homomorphed()
self.symmetry = Symmetry(symmetry) self.symmetry = Symmetry(symmetry)

View File

@ -17,6 +17,7 @@ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options
Add quaternion and/or Bunge Euler angle representation of crystal lattice orientation. Add quaternion and/or Bunge Euler angle representation of crystal lattice orientation.
Orientation is given by quaternion, Euler angles, rotation matrix, or crystal frame coordinates Orientation is given by quaternion, Euler angles, rotation matrix, or crystal frame coordinates
(i.e. component vectors of rotation matrix). (i.e. component vectors of rotation matrix).
Additional (globally fixed) rotations of the lab frame and/or crystal frame can be applied.
""", version = scriptID) """, version = scriptID)
@ -24,24 +25,27 @@ outputChoices = ['quaternion','rodrigues','eulers']
parser.add_option('-o', '--output', parser.add_option('-o', '--output',
dest = 'output', dest = 'output',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'output orientation formats {%s}'%(','.join(outputChoices))) help = 'output orientation formats {{{}}}'.format(', '.join(outputChoices)))
parser.add_option('-r', '--rotation',
dest='rotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis to (pre)rotate orientation')
parser.add_option('-s', '--symmetry', parser.add_option('-s', '--symmetry',
dest = 'symmetry', dest = 'symmetry',
type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string', type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string',
help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'angles are given in degrees [%default]')
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',
dest='crystalrotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis of additional crystal frame rotation')
parser.add_option('-e', '--eulers', parser.add_option('-e', '--eulers',
dest = 'eulers', dest = 'eulers',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'Euler angles label') 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', parser.add_option('-m', '--matrix',
dest = 'matrix', dest = 'matrix',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
@ -65,7 +69,8 @@ parser.add_option('-q', '--quaternion',
parser.set_defaults(output = [], parser.set_defaults(output = [],
symmetry = damask.Symmetry.lattices[-1], symmetry = damask.Symmetry.lattices[-1],
rotation = (0.,1.,1.,1.), # no rotation about 1,1,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, degrees = False,
) )
@ -91,16 +96,16 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
(options.quaternion,4,'quaternion'), (options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested ][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 toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
r = damask.Quaternion().fromAngleAxis(toRadians*options.rotation[0],options.rotation[1:]) # pre-rotation 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 ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name, buffered = False)
buffered = False)
except: continue except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -126,9 +131,9 @@ for name in filenames:
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for output in options.output: for output in options.output:
if output == 'quaternion': table.labels_append(['{}_quat({})'.format( i+1,options.symmetry) for i in xrange(4)]) if output == 'quaternion': table.labels_append(['{}_{}_{}({})'.format(i+1,'quat',options.symmetry,label) for i in xrange(4)])
if output == 'rodrigues': table.labels_append(['{}_rodrigues({})'.format(i+1,options.symmetry) for i in xrange(3)]) elif output == 'rodrigues': table.labels_append(['{}_{}_{}({})'.format(i+1,'rodr',options.symmetry,label) for i in xrange(3)])
if output == 'eulers': table.labels_append(['{}_eulers({})'.format( i+1,options.symmetry) for i in xrange(3)]) elif output == 'eulers': table.labels_append(['{}_{}_{}({})'.format(i+1,'eulr',options.symmetry,label) for i in xrange(3)])
table.head_write() table.head_write()
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------
@ -150,12 +155,12 @@ for name in filenames:
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])), o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])),
symmetry = options.symmetry).reduced() symmetry = options.symmetry).reduced()
o.quaternion = r*o.quaternion o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations
for output in options.output: for output in options.output:
if output == 'quaternion': table.data_append(o.asQuaternion()) if output == 'quaternion': table.data_append(o.asQuaternion())
if output == 'rodrigues': table.data_append(o.asRodrigues()) elif output == 'rodrigues': table.data_append(o.asRodrigues())
if output == 'eulers': table.data_append(o.asEulers('Bunge', degrees=options.degrees)) elif output == 'eulers': table.data_append(o.asEulers('Bunge', degrees=options.degrees))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------