# Calculate velocities from change in deformation of two increments in one step
# Use 'abaqus python report_v.py -h' for help
#
# The progamm assumes [km] and [s] as base units
# and calculates velocities in [mm/a]
# Output:
#   Prints tabular to stdout with the following columns:
#   node, coord_x, coord_y, coord_z, velocity_x, velocity_y, velocity_z
# Return codes:
#   0: No error
#   1: Wrong usage
#   2: Filename error
#   3: ODB database error
# TODO:
#   - Get rid of the required instance name (make this optional)
#   - Define option -o file to redirect output in file
#   - Seperate error messages and data output (use stdin and stderr)
#   - Use standard error codes
#   - Localization
#   - Let the user select other variables then 'U'
#   - Report values instead of velocities also
#
# Version 1.0 (2003-11-07)
# Copyright Kasper D. Fischer <kasper.fischer@ruhr-uni-bochum.de>, 2003

from odbAccess import *
from abaqusConstants import *
import sys, getopt, string, os.path

# Default values
version = '1.0 (2003-11-07)'
authors = 'Kasper D. Fischer <kasper.fischer@ruhr-uni-bochum.de>'
copyyear = '2003'
nsetName = 'SURFACE'
instanceName = None
stepName = None
incNoL = -1
incNoP = -2
odb = None
debug = 0

#Define Usage Function
def usage(text, n):
    if text :
	print text
    print """
Usage:
%s [options] INSTANCE ODBFILE

where
INSTANCE                  the name of the assembly instance to process
                               (mandatory)
ODBFILE                   path and name of the odb-file to process
                               (mandatory)
and options are:
   -n or --nodeset        give the nodeset to process (default: 'SURFACE')
   -f or --firstinc       give the number of the first increment to
                               calculate the velocity
                               (default: second last increment)
   -l or --lastinc        give the number of the last increment to
                               calculate the velocity
                               (default: last increment)
   -s or --step           give the name of the step to use (default: last step)
   -h or --help           show this message
   -d or --debug          print debug information (ommits output of velocities)
   -V or --version        print version and copyright information
         """ %os.path.basename(sys.argv[0])
    sys.exit(n)

# Define Error Function
def error(text,n):
    print text
    if odb:
	odb.close()
    sys.exit(n)

# MAIN
# Get command line arguments.
try:
    optlist, args =getopt.getopt(sys.argv[1:],\
				 'n:f:l:s:hdV', \
				 ['nodeset', 'firstinc', 'lastinc', 'step', 'help', 'debug', 'version'])
    for opt, value in optlist:
	if debug:
	    print 'Processing option %s' % opt
	if opt in ('-n', '--nodeset'):
	    nsetName = value
	elif opt in ('-f', '--firstinc'):
	    incNoP = string.atoi(value)
	elif opt in ('-l', '--lastinc'):
	    incNoL = string.atoi(value)
	elif opt in ('-s', '--step'):
	    stepName = value
	elif opt in ('-h', '--help'):
	    usage('',0)
	elif opt in ('-V', '--version'):
	    print 'This is version %s of %s' % (version, os.path.basename(sys.argv[0]))
	    print 'Copyright %s by %s' % (copyyear, authors)
	    sys.exit(0)
	elif opt in ('-d', '--debug'):
	    debug = 1
	    print 'Processing option %s' % opt
except getopt.GetoptError:
    usage('Invalid option or value!', 1)
except ValueError:
    usage('Invalid option value!', 1)
try:
    instanceName = args[-2]
    infile = args[-1]
except:
    usage('Please specify the odb file and the instance name.', 1)

# open ODB file
if not os.path.isfile(infile):
    error('File %s not found. Check path and name.' %infile, 2)
odb = openOdb(path=infile)
if not odb:
    error('ODB file %s could not be opened.' % infile, 3)
myAssembly = odb.rootAssembly

# Check to see if the instance exists in the assembly and get it.
try:
    myInstance = odb.rootAssembly.instances[instanceName]
except KeyError:
   error('An instance named %s does not exist in the '\
          'output database %s' % (instanceName, infile), 3)

# Check to see if the node set exists in the instance and get it.
try:
   regionNodes = myInstance.nodeSets[nsetName]
except KeyError:
    error('A node set named %s does not exist in the '\
          'instance %s' % (nsetName, instanceName), 3)

# Ceck to see if requested frames of requested (or last) step exist and get it.
try:
    if stepName == None:
        for step in odb.steps.keys():
            stepName = step
    mystep = odb.steps[stepName]
except KeyError:
    error('Step %s does not exists.' % stepName, 3)
try:
    prevFrame = mystep.frames[incNoP]
    lastFrame = mystep.frames[incNoL]
except (KeyError, IndexError):
    error('First (%d) or second increment (%d) not found in step %s.'\
	  % (incNoP, incNoL, stepName), 3)

# get times, deltaTime
try:
    timeL = lastFrame.frameValue
    timeP = prevFrame.frameValue
    deltaTime = (timeL-timeP)
except:
    error('No time information found.',3)
if debug:
    print 'Processing ODB file %s, step %s, increments %d and %d.'\
	  % (infile, stepName, incNoP, incNoL)
    print 'Found increment with time %8.3f' % timeL
    print 'Found increment with time %8.3f' % timeP
    print 'Time difference is %8.3f' % deltaTime

# get data
try:
    dataL = lastFrame.fieldOutputs['U']
    dataP = prevFrame.fieldOutputs['U']
except KeyError:
    error('No data values U found at requested increments', 3)
try:
    regionDataL = dataL.getSubset(region=regionNodes)
    regionDataP = dataP.getSubset(region=regionNodes)
except:
    error('No values found within nodeset %s' %nsetName, 3)

# calculate and print out velocities
try:
    velocity = 3.15e13*(regionDataL - regionDataP)/deltaTime
    for v in velocity.values:
	    node = myInstance.getNodeFromLabel(v.nodeLabel)
	    if not debug:
		print '%d' % v.nodeLabel,
		print ' %8.3f %8.3f %8.3f' % node.coordinate,
		print ' %6.4f %6.4f %6.4f' % v.data
except:
    error('Velocities could not be calculated.', 3)
odb.close()
sys.exit(0)
#END
