DAMASK_EICMD/processing/post/addMises.py

107 lines
3.9 KiB
Python
Raw Normal View History

#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,math,string
import numpy as np
2014-07-25 01:51:18 +05:30
from collections import defaultdict
from optparse import OptionParser
import damask
scriptID = string.replace('$Id$','\n','\\n')
2014-12-19 00:56:52 +05:30
scriptName = os.path.splitext(scriptID.split()[1])[0]
def Mises(what,tensor):
dev = tensor - np.trace(tensor)/3.0*np.eye(3)
symdev = 0.5*(dev+dev.T)
return math.sqrt(np.sum(symdev*symdev.T)*
{
'stress': 3.0/2.0,
'strain': 2.0/3.0,
}[what.lower()])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Add vonMises equivalent values for symmetric part of requested strains and/or stresses.
""", version = scriptID)
parser.add_option('-e','--strain',
dest = 'strain',
action = 'extend', metavar = '<string LIST>',
help = 'heading(s) of columns containing strain tensors')
parser.add_option('-s','--stress',
dest = 'stress',
action = 'extend', metavar = '<string LIST>',
help = 'heading(s) of columns containing stress tensors')
parser.set_defaults(strain = [],
stress = [],
)
(options,filenames) = parser.parse_args()
if len(options.stress+options.strain) == 0:
parser.error('no data column specified...')
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False)
except: continue
table.croak('\033[1m'+scriptName+'\033[0m'+(': '+name if name else ''))
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
items = {
'strain': {'dim': 9, 'shape': [3,3], 'labels':options.strain, 'active':[], 'column': []},
'stress': {'dim': 9, 'shape': [3,3], 'labels':options.stress, 'active':[], 'column': []},
}
errors = []
remarks = []
for type, data in items.iteritems():
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))
table.labels_append('Mises({})'.format(what)) # extend ASCII header with new labels
if remarks != []: table.croak(remarks)
if errors != []:
table.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
2011-12-15 14:23:18 +05:30
table.head_write()
# ------------------------------------------ process data ------------------------------------------
2014-07-25 01:51:18 +05:30
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for type, data in items.iteritems():
for column in data['column']:
table.data_append(Mises(type,
np.array(table.data[column:column+data['dim']],'d').reshape(data['shape'])))
2014-07-25 01:51:18 +05:30
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.close() # close input ASCII table (works for stdin)