223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

replace outdated ASCIItable class

This commit is contained in:
Martin Diehl 2020-05-21 17:16:27 +02:00
parent 7aecc37219
commit 06e4327c0b
1 changed files with 52 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
processing/post/addSchmidfactors.py
View File

@ -2,6 +2,7 @@
import os import os
import sys import sys
from io import StringIO
from optparse import OptionParser from optparse import OptionParser
import numpy as np import numpy as np
@ -102,11 +103,11 @@ Add columns listing Schmid factors (and optional trace vector of selected system
""", version = scriptID) """, version = scriptID)
latticeChoices = ('fcc','bcc','hex') lattice_choices = list(slipSystems.keys())
parser.add_option('-l', parser.add_option('-l',
'--lattice', '--lattice',
dest = 'lattice', type = 'choice', choices = latticeChoices, metavar='string', dest = 'lattice', type = 'choice', choices = lattice_choices, metavar='string',
help = 'type of lattice structure [%default] {}'.format(latticeChoices)) help = 'type of lattice structure [%default] {}'.format(lattice_choices))
parser.add_option('--covera', parser.add_option('--covera',
dest = 'CoverA', type = 'float', metavar = 'float', dest = 'CoverA', type = 'float', metavar = 'float',
help = 'C over A ratio for hexagonal systems [%default]') help = 'C over A ratio for hexagonal systems [%default]')
@ -129,31 +130,29 @@ parser.add_option('-o',
parser.set_defaults(force = (0.0,0.0,1.0), parser.set_defaults(force = (0.0,0.0,1.0),
quaternion='orientation', quaternion='orientation',
normal = None, normal = None,
lattice = latticeChoices[0], lattice = lattice_choices[0],
CoverA = np.sqrt(8./3.), CoverA = np.sqrt(8./3.),
) )
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
if filenames == []: filenames = [None]
force = np.array(options.force) force = np.array(options.force)/np.linalg.norm(options.force)
force /= np.linalg.norm(force)
if options.normal is not None: if options.normal is not None:
normal = np.array(options.normal) normal = np.array(options.normal)/np.linalg.norm(options.ormal)
normal /= np.linalg.norm(normal)
if abs(np.dot(force,normal)) > 1e-3: if abs(np.dot(force,normal)) > 1e-3:
parser.error('stress plane normal not orthogonal to force direction') parser.error('stress plane normal not orthogonal to force direction')
else: else:
normal = force normal = force
slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
slip_normal = np.zeros_like(slip_direction)
if options.lattice in ['bcc','fcc']:
if options.lattice in latticeChoices[:2]:
slip_direction = slipSystems[options.lattice][:,:3] slip_direction = slipSystems[options.lattice][:,:3]
slip_normal = slipSystems[options.lattice][:,3:] slip_normal = slipSystems[options.lattice][:,3:]
elif options.lattice == latticeChoices[2]: elif options.lattice == 'hex':
slip_direction = np.zeros((len(slipSystems['hex']),3),'d')
slip_normal = np.zeros_like(slip_direction)
# convert 4 Miller index notation of hex to orthogonal 3 Miller index notation # convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
for i in range(len(slip_direction)): for i in range(len(slip_direction)):
slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5, slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
@ -165,52 +164,29 @@ elif options.lattice == latticeChoices[2]:
slipSystems['hex'][i,7]/options.CoverA, slipSystems['hex'][i,7]/options.CoverA,
]) ])
slip_direction /= np.tile(np.linalg.norm(slip_direction,axis=1),(3,1)).T slip_direction /= np.linalg.norm(slip_direction,axis=1,keepdims=True)
slip_normal /= np.tile(np.linalg.norm(slip_normal ,axis=1),(3,1)).T slip_normal /= np.linalg.norm(slip_normal, axis=1,keepdims=True)
# --- loop over input files ------------------------------------------------------------------------ labels = ['S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]'
'({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\
if filenames == []: filenames = [None] .format(normal = theNormal, direction = theDirection,
) for theNormal,theDirection in zip(slip_normal,slip_direction)]
for name in filenames: for name in filenames:
try:
table = damask.ASCIItable(name = name)
except IOError:
continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------ table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
table.head_read() o = damask.Rotation.from_quaternion(table.get(options.quaternion))
# ------------------------------------------ sanity checks ---------------------------------------- force = np.broadcast_to(force, o.shape+(3,))
if not table.label_dimension(options.quaternion) == 4: normal = np.broadcast_to(normal,o.shape+(3,))
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion)) slip_direction = np.broadcast_to(slip_direction,o.shape+slip_direction.shape)
table.close(dismiss = True) # close ASCIItable and remove empty file slip_normal = np.broadcast_to(slip_normal, o.shape+slip_normal.shape)
continue S = np.abs(np.einsum('ijk,ik->ij',slip_direction,(o@force))*
np.einsum('ijk,ik->ij',slip_normal, (o@normal)))
column = table.label_index(options.quaternion) for i,label in enumerate(labels):
table.add(label,S[:,i],scriptID+' '+' '.join(sys.argv[1:]))
# ------------------------------------------ assemble header --------------------------------------- table.to_ASCII(sys.stdout if name is None else name)
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append(['S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]'
'({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\
.format(normal = theNormal,
direction = theDirection,
) for theNormal,theDirection in zip(slip_normal,slip_direction)])
table.head_write()
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
o = damask.Rotation(np.array(list(map(float,table.data[column:column+4]))))
table.data_append( np.abs( np.sum(slip_direction * (o * force) ,axis=1) \
* np.sum(slip_normal * (o * normal),axis=1)))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables