DAMASK_EICMD/processing/post/addSchmidfactors.py

#!/usr/bin/env python3

import os
import sys
from optparse import OptionParser

import numpy as np

import damask


scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID   = ' '.join([scriptName,damask.version])

slipSystems = {
'fcc':
    np.array([
    [+0,+1,-1 , +1,+1,+1],
    [-1,+0,+1 , +1,+1,+1],
    [+1,-1,+0 , +1,+1,+1],
    [+0,-1,-1 , -1,-1,+1],
    [+1,+0,+1 , -1,-1,+1],
    [-1,+1,+0 , -1,-1,+1],
    [+0,-1,+1 , +1,-1,-1],
    [-1,+0,-1 , +1,-1,-1],
    [+1,+1,+0 , +1,-1,-1],
    [+0,+1,+1 , -1,+1,-1],
    [+1,+0,-1 , -1,+1,-1],
    [-1,-1,+0 , -1,+1,-1],
   ],'d'),
'bcc':
    np.array([
    [+1,-1,+1 , +0,+1,+1],
    [-1,-1,+1 , +0,+1,+1],
    [+1,+1,+1 , +0,-1,+1],
    [-1,+1,+1 , +0,-1,+1],
    [-1,+1,+1 , +1,+0,+1],
    [-1,-1,+1 , +1,+0,+1],
    [+1,+1,+1 , -1,+0,+1],
    [+1,-1,+1 , -1,+0,+1],
    [-1,+1,+1 , +1,+1,+0],
    [-1,+1,-1 , +1,+1,+0],
    [+1,+1,+1 , -1,+1,+0],
    [+1,+1,-1 , -1,+1,+0],
    [-1,+1,+1 , +2,+1,+1],
    [+1,+1,+1 , -2,+1,+1],
    [+1,+1,-1 , +2,-1,+1],
    [+1,-1,+1 , +2,+1,-1],
    [+1,-1,+1 , +1,+2,+1],
    [+1,+1,-1 , -1,+2,+1],
    [+1,+1,+1 , +1,-2,+1],
    [-1,+1,+1 , +1,+2,-1],
    [+1,+1,-1 , +1,+1,+2],
    [+1,-1,+1 , -1,+1,+2],
    [-1,+1,+1 , +1,-1,+2],
    [+1,+1,+1 , +1,+1,-2],
   ],'d'),
'hex':
    np.array([
    [+2,-1,-1,+0 , +0,+0,+0,+1],
    [-1,+2,-1,+0 , +0,+0,+0,+1],
    [-1,-1,+2,+0 , +0,+0,+0,+1],
    [+2,-1,-1,+0 , +0,+1,-1,+0],
    [-1,+2,-1,+0 , -1,+0,+1,+0],
    [-1,-1,+2,+0 , +1,-1,+0,+0],
    [-1,+1,+0,+0 , +1,+1,-2,+0],
    [+0,-1,+1,+0 , -2,+1,+1,+0],
    [+1,+0,-1,+0 , +1,-2,+1,+0],
    [-1,+2,-1,+0 , +1,+0,-1,+1],
    [-2,+1,+1,+0 , +0,+1,-1,+1],
    [-1,-1,+2,+0 , -1,+1,+0,+1],
    [+1,-2,+1,+0 , -1,+0,+1,+1],
    [+2,-1,-1,+0 , +0,-1,+1,+1],
    [+1,+1,-2,+0 , +1,-1,+0,+1],
    [-2,+1,+1,+3 , +1,+0,-1,+1],
    [-1,-1,+2,+3 , +1,+0,-1,+1],
    [-1,-1,+2,+3 , +0,+1,-1,+1],
    [+1,-2,+1,+3 , +0,+1,-1,+1],
    [+1,-2,+1,+3 , -1,+1,+0,+1],
    [+2,-1,-1,+3 , -1,+1,+0,+1],
    [+2,-1,-1,+3 , -1,+0,+1,+1],
    [+1,+1,-2,+3 , -1,+0,+1,+1],
    [+1,+1,-2,+3 , +0,-1,+1,+1],
    [-1,+2,-1,+3 , +0,-1,+1,+1],
    [-1,+2,-1,+3 , +1,-1,+0,+1],
    [-2,+1,+1,+3 , +1,-1,+0,+1],
    [-1,-1,+2,+3 , +1,+1,-2,+2],
    [+1,-2,+1,+3 , -1,+2,-1,+2],
    [+2,-1,-1,+3 , -2,+1,+1,+2],
    [+1,+1,-2,+3 , -1,-1,+2,+2],
    [-1,+2,-1,+3 , +1,-2,+1,+2],
    [-2,+1,+1,+3 , +2,-1,-1,+2],
   ],'d'),
}

# --------------------------------------------------------------------
#                                MAIN
# --------------------------------------------------------------------

parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
Add columns listing Schmid factors (and optional trace vector of selected system) for given Euler angles.

""", version = scriptID)

latticeChoices = ('fcc','bcc','hex')
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 [%default]')
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',
                  dest = 'normal',
                  type = 'float', nargs = 3, metavar = 'float float float',
                  help = 'stress plane normal in lab frame, per default perpendicular to the force')
parser.add_option('-o',
                  '--orientation',
                  dest = 'quaternion',
                  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 = np.sqrt(8./3.),
                   )

(options, filenames) = parser.parse_args()

force = np.array(options.force)
force /= np.linalg.norm(force)

if options.normal is not None:
  normal = np.array(options.normal)
  normal /= np.linalg.norm(normal)
  if abs(np.dot(force,normal)) > 1e-3:
    parser.error('stress plane normal not orthogonal to force direction')
else:
  normal = force

slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
slip_normal    = np.zeros_like(slip_direction)


if options.lattice in latticeChoices[:2]:
  slip_direction = slipSystems[options.lattice][:,:3]
  slip_normal    = slipSystems[options.lattice][:,3:]
elif options.lattice == latticeChoices[2]:
  # convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
  for i in range(len(slip_direction)):
    slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
                              (slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3),
                               slipSystems['hex'][i,3]*options.CoverA,
                              ])
    slip_normal[i]    = np.array([slipSystems['hex'][i,4],
                              (slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3),
                               slipSystems['hex'][i,7]/options.CoverA,
                              ])

slip_direction /= np.tile(np.linalg.norm(slip_direction,axis=1),(3,1)).T
slip_normal    /= np.tile(np.linalg.norm(slip_normal   ,axis=1),(3,1)).T

# --- loop over input files ------------------------------------------------------------------------

if filenames == []: filenames = [None]

for name in filenames:
  try:
    table = damask.ASCIItable(name = name)
  except IOError:
    continue
  damask.util.report(scriptName,name)

# ------------------------------------------ read header ------------------------------------------

  table.head_read()

# ------------------------------------------ sanity checks ----------------------------------------
  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:]))
  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