From 94efd576639352ec48e723aaecfe7b178a542580 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 16 Aug 2010 20:47:27 +0000 Subject: [PATCH] started to put Marc/Mentat related scripts to corresponding post/pre processing folders. We should move development from my Code folder to here (sorry, lost history then...) --- processing/post/postResults | 725 +++++++++++++++++ processing/pre/marc_addUserOutput | 152 ++++ processing/pre/mentat_colorMap | 180 +++++ .../mentat_patchFromReconstructedBoundaries | 757 ++++++++++++++++++ processing/pre/mentat_pbcOnBoxMesh | 175 ++++ 5 files changed, 1989 insertions(+) create mode 100755 processing/post/postResults create mode 100755 processing/pre/marc_addUserOutput create mode 100755 processing/pre/mentat_colorMap create mode 100755 processing/pre/mentat_patchFromReconstructedBoundaries create mode 100755 processing/pre/mentat_pbcOnBoxMesh diff --git a/processing/post/postResults b/processing/post/postResults new file mode 100755 index 000000000..9cf20c186 --- /dev/null +++ b/processing/post/postResults @@ -0,0 +1,725 @@ +#!/usr/bin/env python + +releases = ['2010b3','2008r1','2007r1','2005r3'] + + +import os, sys, math, re, threading, time +from optparse import OptionParser, OptionGroup, Option, SUPPRESS_HELP +for release in releases: + libPath = '/msc/mentat%s/shlib/'%release + if os.path.exists(libPath): + sys.path.append(libPath) + break +from py_post import * + + +# ----------------------------- +class MyOption(Option): +# ----------------------------- +# used for definition of new option parser action 'extend', which enables to take multiple option arguments +# taken from online tutorial http://docs.python.org/library/optparse.html + + ACTIONS = Option.ACTIONS + ("extend",) + STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",) + TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",) + ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",) + + def take_action(self, action, dest, opt, value, values, parser): + if action == "extend": + lvalue = value.split(",") + values.ensure_value(dest, []).extend(lvalue) + else: + Option.take_action(self, action, dest, opt, value, values, parser) + + +# ----------------------------- +class backgroundMessage(threading.Thread): +# ----------------------------- + + def __init__(self): + threading.Thread.__init__(self) + self.message = '' + self.new_message = '' + self.counter = 0 + self.symbols = ['- ', '\ ', '| ', '/ '] + self.waittime = 0.5 + + def __quit__(self): + length = len(self.message) + len(self.symbols[self.counter]) + sys.stderr.write(chr(8)*length + ' '*length + chr(8)*length) + sys.stderr.write('') + + def run(self): + while not threading.enumerate()[0]._Thread__stopped: + time.sleep(self.waittime) + self.update_message() + self.__quit__() + + def set_message(self, new_message): + self.new_message = new_message + self.print_message() + + def print_message(self): + length = len(self.message) + len(self.symbols[self.counter]) + sys.stderr.write(chr(8)*length + ' '*length + chr(8)*length) # delete former message + sys.stderr.write(self.symbols[self.counter] + self.new_message) # print new message + self.message = self.new_message + + def update_message(self): + self.counter = (self.counter + 1)%len(self.symbols) + self.print_message() + + +# ----------------------------- +def ipCoords(elemType, nodalCoordinates): +# +# returns IP coordinates for a given element +# ----------------------------- + + nodeWeightsPerNode = { + 7: [ [27.0, 9.0, 3.0, 9.0, 9.0, 3.0, 1.0, 3.0], + [ 9.0, 27.0, 9.0, 3.0, 3.0, 9.0, 3.0, 1.0], + [ 3.0, 9.0, 27.0, 9.0, 1.0, 3.0, 9.0, 3.0], + [ 9.0, 3.0, 9.0, 27.0, 3.0, 1.0, 3.0, 9.0], + [ 9.0, 3.0, 1.0, 3.0, 27.0, 9.0, 3.0, 9.0], + [ 3.0, 9.0, 3.0, 1.0, 9.0, 27.0, 9.0, 3.0], + [ 1.0, 3.0, 9.0, 3.0, 3.0, 9.0, 27.0, 9.0], + [ 3.0, 1.0, 3.0, 9.0, 9.0, 3.0, 9.0, 27.0] ], + 117: [ [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], + [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] ], + 136: [ [42.0, 15.0, 15.0, 14.0, 5.0, 5.0], + [15.0, 42.0, 15.0, 5.0, 14.0, 5.0], + [15.0, 15.0, 42.0, 5.0, 5.0, 14.0], + [14.0, 5.0, 5.0, 42.0, 15.0, 15.0], + [ 5.0, 14.0, 5.0, 15.0, 42.0, 15.0], + [ 5.0, 5.0, 14.0, 15.0, 15.0, 42.0] ], + } + + ipCoordinates = [[0.0,0.0,0.0] for i in range(len(nodalCoordinates))] + for ip in range(len(nodeWeightsPerNode[elemType])): + for node in range(len(nodeWeightsPerNode[elemType][ip])): + for i in range(3): + ipCoordinates[ip][i] += nodeWeightsPerNode[elemType][ip][node] * nodalCoordinates[node][i] + for i in range(3): + ipCoordinates[ip][i] /= sum(nodeWeightsPerNode[elemType][ip]) + + return ipCoordinates + + +# ----------------------------- +def sortBySeparation(dataArray, criteria, offset): +# +# sorting of groupValue array according to list of criteria +# ----------------------------- + where = { + 'elem': 1, + 'node': 2, + 'grain': 3, + 'x': 4, + 'y': 5, + 'z': 6, + } + + theKeys = [] + for criterium in criteria: + if criterium in where: + theKeys.append('x[%i]'%(offset+where[criterium])) + exec('sortedArray = sorted(dataArray,key=lambda x:(%s))'%(','.join(theKeys))) + return sortedArray + + +# ----------------------------- +def substituteLocation(string, mesh, coords): +# +# do variable interpolation in group and filter strings +# ----------------------------- + substitute = string + substitute = substitute.replace('elem', str(mesh[0])) + substitute = substitute.replace('node', str(mesh[1])) + substitute = substitute.replace('grain', str(mesh[2])) + substitute = substitute.replace('x', '%.6g'%coords[0]) + substitute = substitute.replace('y', '%.6g'%coords[1]) + substitute = substitute.replace('z', '%.6g'%coords[2]) + return substitute + + +# ----------------------------- +def average(theList): +# +# calcs the average of a list of numbers +# ----------------------------- + + return sum(map(float,theList))/len(theList) + + +# ----------------------------- +def mapFunc(label, chunks, func): +# +# applies the function defined by "func" +# (can be either 'min','max','avg', 'sum', or user specified) +# to a list of lists of data +# ----------------------------- + + illegal = { + 'eulerangles': ['min','max','avg','sum'], + 'defgrad': ['min','max','avg','sum'], + 'orientation': ['min','max','sum'], + } + + if label.lower() in illegal and func in illegal[label.lower()]: # for illegal mappings:... + return ['n/a' for i in range(len(chunks[0]))] # ...return 'n/a' + + else: + if func in ['min','max','avg']: + mapped = [{ 'min': lambda x: min(x), + 'max': lambda x: max(x), + 'avg': lambda x: average(x), + 'sum': lambda x: sum(x), + }[func](column) for column in zip(*chunks)] # map one of the standard functions to colums in chunks + if label.lower() == 'orientation': # orientation is special case:... + orientationNorm = math.sqrt(sum([q*q for q in mapped])) # ...calc norm of average quaternion + mapped = map(lambda x: x/orientationNorm, mapped) # ...renormalize quaternion + else: + try: + mapped = eval('map(%s,zip(*chunks))'%func) # map user defined function to colums in chunks + except: + mapped = ['n/a' for i in range(len(chunks[0]))] + + return mapped + + +# ----------------------------- +def OpenPostfile(name): +# +# open postfile with extrapolation mode "translate" +# ----------------------------- + + p = post_open(name) + p.extrapolation('translate') + p.moveto(1) + + return p + + +# ----------------------------- +def ParseOutputFormat(filename,what,me): +# +# parse .output* files in order to get a list of outputs +# ----------------------------- + + format = {'outputs':{},'specials':{'brothers':[]}} + for prefix in ['']+map(str,range(1,17)): + if os.path.exists(prefix+filename+'.output'+what): + break + try: + file = open(prefix+filename+'.output'+what) + content = file.readlines() + file.close() + except: + return format + + tag = '' + tagID = 0 + for line in content: + if re.match("\s*$",line) or re.match("#",line): # skip blank lines and comments + continue + m = re.match("\[(.+)\]",line) # look for block indicator + if m: # next section + tag = m.group(1) + tagID += 1 + format['specials']['brothers'].append(tag) + if tag == me or (me.isdigit() and tagID == int(me)): + format['specials']['_id'] = tagID + format['outputs'] = [] + tag = me + else: # data from section + if tag == me: + (output,length) = line.split() + output.lower() + if length.isdigit(): + length = int(length) + if re.match("\((.+)\)",output): # special data, e.g. (Ngrains) + format['specials'][output] = length + elif length > 0: + format['outputs'].append([output,length]) + return format + + +# ----------------------------- +def ParsePostfile(p,filename, outputFormat): +# +# parse postfile in order to get position and labels of outputs +# needs "outputFormat" for mapping of output names to postfile output indices +# ----------------------------- + + # --- build statistics + + stat = { \ + 'IndexOfLabel': {}, \ + 'Title': p.title(), \ + 'Extrapolation': p.extrapolate, \ + 'NumberOfIncrements': p.increments(), \ + 'NumberOfNodes': p.nodes(), \ + 'NumberOfNodalScalars': p.node_scalars(), \ + 'LabelOfNodalScalar': [None]*p.node_scalars() , \ + 'NumberOfElements': p.elements(), \ + 'NumberOfElementalScalars': p.element_scalars(), \ + 'LabelOfElementalScalar': [None]*p.element_scalars() , \ + 'NumberOfElementalTensors': p.element_tensors(), \ + 'LabelOfElementalTensor': [None]*p.element_tensors(), \ + } + + # --- find labels + + for labelIndex in range(stat['NumberOfNodalScalars']): + label = p.node_scalar_label(labelIndex) + stat['IndexOfLabel'][label] = labelIndex + stat['LabelOfNodalScalar'][labelIndex] = label + + for labelIndex in range(stat['NumberOfElementalScalars']): + label = p.element_scalar_label(labelIndex) + stat['IndexOfLabel'][label] = labelIndex + stat['LabelOfElementalScalar'][labelIndex] = label + + for labelIndex in range(stat['NumberOfElementalTensors']): + label = p.element_tensor_label(labelIndex) + stat['IndexOfLabel'][label] = labelIndex + stat['LabelOfElementalTensor'][labelIndex] = label + + if 'User Defined Variable 1' in stat['IndexOfLabel']: + stat['IndexOfLabel']['GrainCount'] = stat['IndexOfLabel']['User Defined Variable 1'] + + if 'GrainCount' in stat['IndexOfLabel']: # does the result file contain relevant user defined output at all? + startIndex = stat['IndexOfLabel']['GrainCount'] - 1 + + # We now have to find a mapping for each output label as defined in the .output* files to the output position in the post file + # Since we know where the user defined outputs start ("startIndex"), we can simply assign increasing indices to the labels + # given in the .output* file + + offset = 0 + stat['LabelOfElementalScalar'][startIndex + 2 + offset] = 'HomogenizationCount' + for var in outputFormat['Homogenization']['outputs']: + if var[1] > 1: + for i in range(var[1]): + stat['IndexOfLabel']['%i_%s'%(i+1,var[0])] = startIndex + 2 + offset + (i+1) + else: + stat['IndexOfLabel']['%s'%(var[0])] = startIndex + 2 + offset + 1 + offset += var[1] + + for grain in range(outputFormat['Homogenization']['specials']['(ngrains)']): + stat['IndexOfLabel']['%i_CrystalliteCount'%(grain+1)] = startIndex + 3 + offset + for var in outputFormat['Crystallite']['outputs']: + if var[1] > 1: + for i in range(var[1]): + stat['IndexOfLabel']['%i_%i_%s'%(grain+1,i+1,var[0])] = startIndex + 3 + offset + (i+1) + else: + stat['IndexOfLabel']['%i_%s'%(grain+1,var[0])] = startIndex + 3 + offset + 1 + offset += var[1] + + stat['IndexOfLabel']['%i_ConstitutiveCount'%(grain+1)] = startIndex + 4 + offset + for var in outputFormat['Constitutive']['outputs']: + if var[1] > 1: + for i in range(var[1]): + stat['IndexOfLabel']['%i_%i_%s'%(grain+1,i+1,var[0])] = startIndex + 4 + offset + (i+1) + else: + stat['IndexOfLabel']['%i_%s'%(grain+1,var[0])] = startIndex + 4 + offset + 1 + offset += var[1] + + return stat + + +# ----------------------------- +def SummarizePostfile(stat,where=sys.stdout): +# ----------------------------- + + where.write('title:\t%s'%stat['Title'] + '\n\n') + where.write('extraplation:\t%s'%stat['Extrapolation'] + '\n\n') + where.write('increments:\t%i+1'%(stat['NumberOfIncrements']-1) + '\n\n') + where.write('nodes:\t%i'%stat['NumberOfNodes'] + '\n\n') + where.write('elements:\t%i'%stat['NumberOfElements'] + '\n\n') + where.write('nodal scalars:\t%i'%stat['NumberOfNodalScalars'] + '\n\n ' + '\n '.join(stat['LabelOfNodalScalar']) + '\n\n') + where.write('elemental scalars:\t%i'%stat['NumberOfElementalScalars'] + '\n\n ' + '\n '.join(stat['LabelOfElementalScalar']) + '\n\n') + where.write('elemental tensors:\t%i'%stat['NumberOfElementalTensors'] + '\n\n ' + '\n '.join(stat['LabelOfElementalTensor']) + '\n\n') + + return True + + +# ----------------------------- +# MAIN FUNCTION STARTS HERE +# ----------------------------- + +# --- input parsing + +parser = OptionParser(option_class=MyOption, usage='%prog [options] resultfile', description = """ +Extract data from a .t16 MSC.Marc results file. + +List of output variables is given by options '--ns','--es','--et','--ho','--cr','--co'. + +Filter and separations use 'elem','node','grain', and 'x','y','z' as key words. +Example: +1) get averaged results in slices perpendicular to x for all positive y coordinates +--filter 'y >= 0.0' --separation x --map 'avg' +2) global sum of squared data falling into first quadrant arc between R1 and R2 +--filter 'x*x + y*y >= R1*R1 and x*x + y*y <= R2*R2' --map 'lambda list: sum([item*item for item in list])' + +$Id: postResults 205 2010-06-08 15:23:31Z MPIE\p.eisenlohr $ +""") + +parser.add_option('-i','--info', action='store_true', dest='info', \ + help='list contents of resultfile [%default]') +parser.add_option('-d','--dir', dest='directory', \ + help='name of subdirectory to hold output [%default]') +parser.add_option('-r','--range', dest='range', type='int', nargs=3, \ + help='range of increments to output (start, end, step) [all]') +parser.add_option('-m','--map', dest='func', type='string', \ + help='data reduction mapping ["%default"] out of min, max, avg, sum or user-lambda') + +group_material = OptionGroup(parser,'Material identifier') +group_special = OptionGroup(parser,'Special outputs') +group_general = OptionGroup(parser,'General outputs') + +group_material.add_option('--homogenization', dest='homog', type='string', \ + help='homogenization identifier (as string or integer [%default])') +group_material.add_option('--crystallite', dest='cryst', type='string', \ + help='crystallite identifier (as string or integer [%default])') +group_material.add_option('--phase', dest='phase', type='string', \ + help='phase identifier (as string or integer [%default])') + +group_special.add_option('-t','--time', action='store_true', dest='time', \ + help='output time of increment [%default]') +group_special.add_option('-f','--filter', dest='filter', type='string', \ + help='condition(s) to filter results [%default]') +group_special.add_option('--separation', action='extend', dest='separation', type='string', \ + help='properties to separate results [%default]') +parser.add_option('-s','--split', action='store_true', dest='separateFiles', \ + help='split output per increment [%default]') + +group_general.add_option('--ns', action='extend', dest='nodalScalar', type='string', \ + help='list of nodal scalars to extract') +group_general.add_option('--es', action='extend', dest='elementalScalar', type='string', \ + help='list of elemental scalars to extract') +group_general.add_option('--et', action='extend', dest='elementalTensor', type='string', \ + help='list of elemental tensors to extract') +group_general.add_option('--ho', action='extend', dest='homogenizationResult', type='string', \ + help='list of homogenization results to extract') +group_general.add_option('--cr', action='extend', dest='crystalliteResult', type='string', \ + help='list of crystallite results to extract') +group_general.add_option('--co', action='extend', dest='constitutiveResult', type='string', \ + help='list of constitutive results to extract') + +parser.add_option_group(group_material) +parser.add_option_group(group_general) +parser.add_option_group(group_special) + +parser.set_defaults(info = False) +parser.set_defaults(directory = 'postProc') +parser.set_defaults(func = 'avg') +parser.set_defaults(homog = '1') +parser.set_defaults(cryst = '1') +parser.set_defaults(phase = '1') +parser.set_defaults(filter = '') +parser.set_defaults(separation = []) +parser.set_defaults(inc = False) +parser.set_defaults(time = False) +parser.set_defaults(separateFiles = False) + +(options, file) = parser.parse_args() + +bg = backgroundMessage() +bg.start() + + +# --- sanity checks + +if not file: + parser.print_help() + parser.error('no file specified...') + +if options.constitutiveResult and not options.phase: + parser.print_help() + parser.error('constitutive results require phase...') + +if options.nodalScalar and ( options.elementalScalar or options.elementalTensor + or options.homogenizationResult or options.crystalliteResult or options.constitutiveResult ): + parser.print_help() + parser.error('not allowed to mix nodal with elemental results...') + + +# --- parse .output and .t16 files + +bg.set_message('parsing .output and .t16 files...') + +filename = os.path.splitext(file[0])[0] +dirname = os.path.abspath(os.path.dirname(filename))+os.sep+options.directory +if not os.path.isdir(dirname): + os.mkdir(dirname,0755) + +outputFormat = {} +me = { + 'Homogenization': options.homog, + 'Crystallite': options.cryst, + 'Constitutive': options.phase, + } +for what in me: + outputFormat[what] = ParseOutputFormat(filename, what, me[what]) + if not '_id' in outputFormat[what]['specials']: + print "'%s' not found in <%s>"%(me[what], what) + print '\n'.join(map(lambda x:' '+x, outputFormat[what]['specials']['brothers'])) + sys.exit(1) + +p = OpenPostfile(filename+'.t16') +stat = ParsePostfile(p, filename, outputFormat) + + +# --- sanity check for output variables +# for mentat variables (nodalScalar,elementalScalar,elementalTensor) we simply have to check whether the label is found in the stat[indexOfLabel] dictionary +# for user defined variables (homogenizationResult,crystalliteResult,constitutiveResult) we have to check the corresponding outputFormat, since the namescheme in stat['IndexOfLabel'] is different + +for opt in ['nodalScalar','elementalScalar','elementalTensor','homogenizationResult','crystalliteResult','constitutiveResult']: + if eval('options.%s'%opt): + for label in eval('options.%s'%opt): + if (opt in ['nodalScalar','elementalScalar','elementalTensor'] and not label in stat['IndexOfLabel']) \ + or (opt in ['homogenizationResult','crystalliteResult','constitutiveResult'] \ + and (not outputFormat[opt[:-6].capitalize()]['outputs'] or not label in zip(*outputFormat[opt[:-6].capitalize()]['outputs'])[0])): + parser.error('%s "%s" unknown...'%(opt,label)) + + +# --- output info + +if options.info: + print '\nMentat release %s\n'%release + SummarizePostfile(stat,sys.stderr) + + print '\nUser Defined Outputs' + for what in me: + print '\n ',what,':' + for output in outputFormat[what]['outputs']: + print ' ',output + + sys.exit(0) + + +# --- get output data from .t16 file + + +if options.range: + increments = range( max(0,options.range[0]), + min(stat['NumberOfIncrements'],options.range[1]+1), + options.range[2]) +else: + increments = range(stat['NumberOfIncrements']-1) + +fileOpen = False +assembleHeader = True +header = [] + +for increment in increments: + p.moveto(increment+1) + bg.set_message('read data from increment %i...'%increment) + data = {} + + if options.nodalScalar: + for n in range(stat['NumberOfNodes']): + nodeID = p.node_id(n) + nodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z] + elemID = 0 + grainID = 0 + + # --- filter valid locations + + filter = substituteLocation(options.filter, [elemID,nodeID,grainID], nodeCoordinates) # generates an expression that is only true for the locations specified by options.filter + if filter != '' and not eval(filter): # for all filter expressions that are not true:... + continue # ... ignore this data point and continue with next + + # --- group data locations + + group = substituteLocation('#'.join(options.separation), [elemID,nodeID,grainID], nodeCoordinates) # generates a unique key for a group of separated data based on the separation criterium for the location + if group not in data: # create a new group if not yet present + data[group] = [] + data[group].append([]) # append a new list for each group member; each list will contain dictionaries with keys 'label, and 'content' for the associated data + data[group][-1].append({ + 'label': 'location', + 'content': [elemID,nodeID,grainID] + nodeCoordinates, + }) # first entry in this list always contains the location data + + # --- get data from t16 file + + for label in options.nodalScalar: + if assembleHeader: + header.append(label.replace(' ','')) + data[group][-1].append({ + 'label': label, + 'content': [ p.node_scalar(n,stat['IndexOfLabel'][label]) ], + }) + + assembleHeader = False + + else: + for e in range(stat['NumberOfElements']): + nodeCoordinates = map(lambda node: [node.x, node.y, node.z], map(p.node, map(p.node_sequence,p.element(e).items))) + ipCoordinates = ipCoords(p.element(e).type, nodeCoordinates) + elemID = p.element_id(e) + for n in range(p.element(e).len): + nodeID = p.element(e).items[n] + for g in range(('GrainCount' in stat['IndexOfLabel'] and int(p.element_scalar(e, stat['IndexOfLabel']['GrainCount'])[0].value)) + or 1): + grainID = g + 1 + + # --- filter valid locations + + filter = substituteLocation(options.filter, [elemID,nodeID,grainID], ipCoordinates[n]) # generates an expression that is only true for the locations specified by options.filter + if filter != '' and not eval(filter): # for all filter expressions that are not true:... + continue # ... ignore this data point and continue with next + + # --- group data locations + + group = substituteLocation('#'.join(options.separation), [elemID,nodeID,grainID], ipCoordinates[n]) # generates a unique key for a group of separated data based on the separation criterium for the location + if group not in data: # create a new group if not yet present + data[group] = [] + data[group].append([]) # append a new list for each group member; each list will contain dictionaries with keys 'label, and 'content' for the associated data + data[group][-1].append({ + 'label': 'location', + 'content': [elemID,nodeID,grainID] + ipCoordinates[n], + }) # first entry in this list always contains the location data + + # --- get data from t16 file + + if options.elementalScalar: + for label in options.elementalScalar: + if assembleHeader: + header.append(label.replace(' ','')) + data[group][-1].append({ + 'label': label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel'][label])[n].value ], + }) + + if options.elementalTensor: + for label in options.elementalTensor: + if assembleHeader: + header += ['%s.%s'%(label.replace(' ',''),component) for component in ['intensity','t11','t22','t33','t12','t23','t13']] + data[group][-1].append({ + 'label': label, + 'content': [ eval("p.element_tensor(e,stat['IndexOfLabel'][label])[n].%s"%component) + for component in ['intensity','t11','t22','t33','t12','t23','t13'] ], + }) + + if options.homogenizationResult: + for label in options.homogenizationResult: + outputIndex = list(zip(*outputFormat['Homogenization']['outputs'])[0]).index(label) # find the position of this output in the outputFormat + length = int(outputFormat['Homogenization']['outputs'][outputIndex][1]) + if length > 1: + if assembleHeader: + header += ['%i_%s'%(component+1,label) for component in range(length)] + data[group][-1].append({ + 'label': label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel']['%i_%s'%(component+1,label)])[n].value + for component in range(length) ], + }) + else: + if assembleHeader: + header.append(label) + data[group][-1].append({ + 'label': label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel']['%s'%label])[n].value ], + }) + + if options.crystalliteResult: + for label in options.crystalliteResult: + outputIndex = list(zip(*outputFormat['Crystallite']['outputs'])[0]).index(label) # find the position of this output in the outputFormat + length = int(outputFormat['Crystallite']['outputs'][outputIndex][1]) + if length > 1: + if assembleHeader: + header += ['%i_%i_%s'%(g+1,component+1,label) for component in range(length)] + data[group][-1].append({ + 'label': label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel']['%i_%i_%s'%(g+1,component+1,label)])[n].value + for component in range(length) ], + }) + else: + if assembleHeader: + header.append('%i_%s'%(g+1,label)) + data[group][-1].append({ + 'label':label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel']['%i_%s'%(g+1,label)])[n].value ], + }) + + if options.constitutiveResult: + for label in options.constitutiveResult: + outputIndex = list(zip(*outputFormat['Constitutive']['outputs'])[0]).index(label) # find the position of this output in the outputFormat + length = int(outputFormat['Constitutive']['outputs'][outputIndex][1]) + if length > 1: + if assembleHeader: + header += ['%i_%i_%s'%(g+1,component+1,label) for component in range(length)] + data[group][-1].append({ + 'label':label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel']['%i_%i_%s'%(g+1,component+1,label)])[n].value + for component in range(length) ], + }) + else: + if assembleHeader: + header.append('%i_%s'%(g+1,label)) + data[group][-1].append({ + 'label':label, + 'content': [ p.element_scalar(e,stat['IndexOfLabel']['%i_%s'%(g+1,label)])[n].value ], + }) + + assembleHeader = False + + if options.separateFiles: + if fileOpen: + file.close() + fileOpen = False + outFilename = eval('"'+eval("'%%s_inc%%0%ii.txt'%(math.log10(max(increments))+1)")+'"%(dirname + os.sep + os.path.split(filename)[1],increment)') + else: + outFilename = '%s.txt'%(dirname + os.sep + os.path.split(filename)[1]) + + # --- write header to file + + if not fileOpen: + file = open(outFilename,'w') + fileOpen = True + file.write('2\theader\n') + file.write('$Id: postResults 205 2010-06-08 15:23:31Z MPIE\p.eisenlohr $\n') + if options.time: + basic = ['inc','time'] + else: + basic = ['inc'] + if options.nodalScalar: + file.write('\t'.join(basic + ['elem','node','grain','node.x','node.y','node.z'] + header) + '\n') + else: + file.write('\t'.join(basic + ['elem','node','grain','ip.x','ip.y','ip.z'] + header) + '\n') + + # --- write data to file + + output = [] + for group in data: + if options.time: + output.append([increment, p.time]) + else: + output.append([increment]) + for chunk in range(len(data[group][0])): + label = data[group][0][chunk]['label'] # name of chunk (e.g. 'orientation', or 'flow stress') + groupContent = [data[group][member][chunk]['content'] for member in range(len(data[group]))] # list of each member's chunk + if label == 'location': + condensedGroupContent = mapFunc(label, groupContent, 'avg') # always average location + if len(groupContent) > 1: # e,n,g nonsense if averaged over more than one entry... + condensedGroupContent[:3] = ['n/a']*3 # ...so return 'n/a' + else: + condensedGroupContent = mapFunc(label, groupContent, options.func) # map function to groupContent to get condensed data of this group's chunk + output[-1] += condensedGroupContent + + for groupvalues in sortBySeparation(output, options.separation, int(options.time)): # sort output according to separation criteria + file.write('\t'.join(map(str,groupvalues)) + '\n') + +if fileOpen: + file.close() + + +# --------------------------- DONE -------------------------------- + diff --git a/processing/pre/marc_addUserOutput b/processing/pre/marc_addUserOutput new file mode 100755 index 000000000..e418e2d9c --- /dev/null +++ b/processing/pre/marc_addUserOutput @@ -0,0 +1,152 @@ +#!/usr/bin/env python +''' +Writes meaningful labels to the marc input file (*.dat) +based on the files +.output +that are written during the first run of the model. +''' +import sys,os,re +from optparse import OptionParser + +# ----------------------------- +def ParseOutputFormat(filename,what,me): +# ----------------------------- + format = {'outputs':{},'specials':{'brothers':[]}} + + outputmetafile = filename+'.output'+what + try: + file = open(outputmetafile) + except: + print('Could not open file %s'%outputmetafile) + raise + else: + content = file.readlines() + file.close() + + tag = '' + tagID = 0 + for line in content: + if re.match("\s*$",line) or re.match("#",line): # skip blank lines and comments + continue + m = re.match("\[(.+)\]",line) # look for block indicator + if m: # next section + tag = m.group(1) + tagID += 1 + format['specials']['brothers'].append(tag) + if tag == me or (me.isdigit() and tagID == int(me)): + format['specials']['_id'] = tagID + format['outputs'] = [] + tag = me + else: # data from section + if tag == me: + (output,length) = line.split() + output.lower() + if length.isdigit(): + length = int(length) + if re.match("\((.+)\)",output): # special data, (e.g. (Ngrains) + format['specials'][output] = length + elif length > 0: + format['outputs'].append([output,length]) + return format + + +parser = OptionParser(usage='%prog [options] Marc.inputfile(s)', description=""" +what I do... +""") +parser.add_option('-n','--number', dest='number', type='int', \ + help='maximum requested User Defined Variable [%default]') +parser.add_option('--homogenization', dest='homog', type='string', \ + help='homogenization identifier (as string or integer [%default])') +parser.add_option('--crystallite', dest='cryst', type='string', \ + help='crystallite identifier (as string or integer [%default])') +parser.add_option('--phase', dest='phase', type='string', \ + help='phase identifier (as string or integer [%default])') +parser.add_option('--use', dest='useFile', type='string', \ + help='Optionally parse output descriptors from '+ + 'different .outputZZZ file. Saves the effort '+ + 'to start a calculation for each job [%default])') +parser.set_defaults(number = 0) +parser.set_defaults(homog = '1') +parser.set_defaults(cryst = '1') +parser.set_defaults(phase = '1') +parser.set_defaults(useFile = '') + +(options, files) = parser.parse_args() + +if not files: + parser.print_help() + parser.error('no file(s) specified...') + +me = { 'Homogenization': options.homog, + 'Crystallite': options.cryst, + 'Constitutive': options.phase, + } + + +for file in files: + if options.useFile != '': + formatFile = os.path.splitext(options.useFile)[0] + else: + formatFile = os.path.splitext(file)[0] + file = os.path.splitext(file)[0]+'.dat' + if not os.path.lexists(file): + print file,'not found' + continue + + print('Scanning format files of: %s'%formatFile) + + if options.number < 1: + outputFormat = {} + + for what in me: + outputFormat[what] = ParseOutputFormat(formatFile,what,me[what]) + if not '_id' in outputFormat[what]['specials']: + print "'%s' not found in <%s>"%(me[what],what) + print '\n'.join(map(lambda x:' '+x,outputFormat[what]['specials']['brothers'])) + sys.exit(1) + + UserVars = ['GrainCount'] + + UserVars += ['HomogenizationCount'] + for var in outputFormat['Homogenization']['outputs']: + if var[1] > 1: + UserVars += ['%i_%s'%(i+1,var[0]) for i in range(var[1])] + else: + UserVars += ['%s'%(var[0]) for i in range(var[1])] + for grain in range(outputFormat['Homogenization']['specials']['(ngrains)']): + UserVars += ['%i_CrystalliteCount'%(grain+1)] + for var in outputFormat['Crystallite']['outputs']: + if var[1] > 1: + UserVars += ['%i_%i_%s'%(grain+1,i+1,var[0]) for i in range(var[1])] + else: + UserVars += ['%i_%s'%(grain+1,var[0]) for i in range(var[1])] + + UserVars += ['%i_ConstitutiveCount'%(grain+1)] + for var in outputFormat['Constitutive']['outputs']: + if var[1] > 1: + UserVars += ['%i_%i_%s'%(grain+1,i+1,var[0]) for i in range(var[1])] + else: + UserVars += ['%i_%s'%(grain+1,var[0]) for i in range(var[1])] + +# Now change *.dat file(s) + print('Adding labels to: %s'%file) + inFile = open(file) + input = inFile.readlines() + inFile.close() + output = open(file,'w') + thisSection = '' + for line in input: + m = re.match('(\w+)\s',line) + if m: + lastSection = thisSection + thisSection = m.group(1) + if (lastSection == 'post' and thisSection == 'parameters'): + if options.number > 0: + for i in range(options.number): + output.write('%10i%10i\n'%(-i-1,0)) + else: + for i in range(len(UserVars)): + output.write('%10i%10i%s\n'%(-i-1,0,UserVars[i])) + if (thisSection != 'post' or not re.match('\s*\-',line)): + output.write(line) + output.close() diff --git a/processing/pre/mentat_colorMap b/processing/pre/mentat_colorMap new file mode 100755 index 000000000..704676f9c --- /dev/null +++ b/processing/pre/mentat_colorMap @@ -0,0 +1,180 @@ +#!/usr/bin/env python + + +import sys, os +from colorsys import * +from optparse import OptionParser + +releases = ['2010','2008r1','2007r1','2005r3'] + +for release in releases: + libPath = '/msc/mentat%s/shlib/'%release + if os.path.exists(libPath): + sys.path.append(libPath) + for subdir in [os.path.join(libPath,file) + for file in os.listdir(libPath) + if os.path.isdir(os.path.join(libPath,file))]: + sys.path.append(subdir) + break +from py_mentat import * + + + +# ----------------------------- +def outMentat(cmd,locals): + if cmd[0:3] == '(!)': + exec(cmd[3:]) + elif cmd[0:3] == '(?)': + cmd = eval(cmd[3:]) + py_send(cmd) + else: + py_send(cmd) + return + + + +# ----------------------------- +def outStdout(cmd,locals): + if cmd[0:3] == '(!)': + exec(cmd[3:]) + elif cmd[0:3] == '(?)': + cmd = eval(cmd[3:]) + print cmd + else: + print cmd + return + + + +# ----------------------------- +def output(cmds,locals,dest): + for cmd in cmds: + if isinstance(cmd,list): + output(cmd,locals,dest) + else: + {\ + 'Mentat': outMentat,\ + 'Stdout': outStdout,\ + }[dest](cmd,locals) + return + + + +# ----------------------------- +def lever(val0, val1, x): + return val0 + (val1 - val0) * x + + + +# ----------------------------- +def symlever(comp, val0, val1, x): + if comp == "hue": + return lever(val0, val1, x) + + if comp == "lightness": + val_middle = max(0.9, val0, val1) + elif comp == "saturation": + val_middle = min(0.1, val0, val1) + if x < 0.5: + return lever(val0, val_middle, 2*x) + else: + return lever(val_middle, val1, 2*x-1) + + + +# ----------------------------- +def colorMap(colors): + cmds = [ "*color %i %f %f %f"%(idx+32,color[0],color[1],color[2]) + for idx,color in enumerate(colors) ] + return cmds + + +# ----------------------------- +# MAIN FUNCTION STARTS HERE +# ----------------------------- + +parser = OptionParser(usage="%prog [options] lower_hls upper_hls", description = """ +Changes the color map in mentat. + +Interpolates colors between "lower_hls" and "upper_hls". +For symmetric scales use option "-s". + +Example colors: +- Non-symmetric scales: 0.167,0.9,0.1 0.167,0.1,0.9 +- Symmetric scales: 0,0.2,0.9 0.333,0.2,0.9 +""") + +parser.add_option("-s","--symmetric", action = "store_true", + dest = "symmetric", \ + help = "symmetric legend [%default]") +parser.add_option("-p", "--port", type = "int",\ + dest = "port",\ + help = "Mentat connection port [%default]") +parser.add_option("-v", "--verbose", action="store_true",\ + dest = "verbose",\ + help = "write Mentat command stream also to stdout [%default]") +parser.set_defaults(symmetric = False) +parser.set_defaults(port=40007) +parser.set_defaults(verbose=False) + + +(options, vars) = parser.parse_args() + + + +### read hlsColors and check if they are valid hls values + +try: + hlsColor_bounds = [[],[]] + for i in range(2): + hlsColor_bounds[i] = map(float, vars[i].split(",")) + + if len(hlsColor_bounds[i]) <> 3: + raise + + hlsColors_limits = [[0,0,0],[1,1,1]] + for j in range(3): + if hlsColor_bounds[i][j] < hlsColors_limits[0][j] or hlsColor_bounds[i][j] > hlsColors_limits[1][j]: + raise + +except: + parser.error("give lower and upper hlsColor as comma separated values") + + + +### interpolate hls values + +nColors = 32 +if options.symmetric: + hlsColors = [ [ symlever(comp, hlsColor_bounds[0][j], hlsColor_bounds[1][j], float(idx)/(nColors-1)) + for j,comp in enumerate(["hue","lightness","saturation"]) ] + for idx in range(nColors) ] +else: + hlsColors = [ [ lever(hlsColor_bounds[0][j], hlsColor_bounds[1][j], float(idx)/(nColors-1)) + for j,comp in enumerate(["hue","lightness","saturation"]) ] + for idx in range(nColors) ] + + + +### convert to rgb values + +rgbColors = [ hls_to_rgb(hlsColor[0], hlsColor[1], hlsColor[2]) + for hlsColor in hlsColors ] + + + +### connect to mentat and change colorMap + +outputLocals = {} +print 'waiting to connect...' +py_connect('',options.port) +print 'connected...' + +cmds = colorMap(rgbColors) +output(cmds,outputLocals,'Mentat') +py_disconnect() + +if options.verbose: + output(cmds,outputLocals,'Stdout') + + \ No newline at end of file diff --git a/processing/pre/mentat_patchFromReconstructedBoundaries b/processing/pre/mentat_patchFromReconstructedBoundaries new file mode 100755 index 000000000..a5aba98d2 --- /dev/null +++ b/processing/pre/mentat_patchFromReconstructedBoundaries @@ -0,0 +1,757 @@ +#!/usr/bin/env python + +releases = ['2010r1','2008r1','2007r1','2005r3'] + + +import sys,os,pwd,math,re +#import Image,ImageDraw +from optparse import OptionParser + +for release in releases: + libPath = '/msc/mentat%s/shlib/'%release + if os.path.exists(libPath): + sys.path.append(libPath) + break +from py_mentat import * + + +def outMentat(cmd,locals): + if cmd[0:3] == '(!)': + exec(cmd[3:]) + elif cmd[0:3] == '(?)': + cmd = eval(cmd[3:]) + py_send(cmd) + else: + py_send(cmd) + return + +def outStdout(cmd,locals): + if cmd[0:3] == '(!)': + exec(cmd[3:]) + elif cmd[0:3] == '(?)': + cmd = eval(cmd[3:]) + print cmd + else: + print cmd + return + + +def output(cmds,locals,dest): + for cmd in cmds: + if isinstance(cmd,list): + output(cmd,locals,dest) + else: + {\ + 'Mentat': outMentat,\ + 'Stdout': outStdout,\ + }[dest](cmd,locals) + return + + +def rcbOrientationParser(content): + + grains = [] + myOrientation = [0.0,0.0,0.0] + for line in content: + if line[0] != '#': # skip comments + for grain in range(2): + myID = int(line.split()[12+grain]) # get grain id + myOrientation = map(float,line.split())[3*grain:3+3*grain] # get orientation + if len(grains) < myID: + for i in range(myID-len(grains)): # extend list to necessary length + grains.append([0.0,0.0,0.0]) + grains[myID-1] = myOrientation # store Euler angles + + return grains + +def rcbParser(content,size,tolerance,imagename,imagesize,border): # parser for TSL-OIM reconstructed boundary files + +# find bounding box + + boxX = [1.*sys.maxint,-1.*sys.maxint] + boxY = [1.*sys.maxint,-1.*sys.maxint] + x = [0.,0.] + y = [0.,0.] + for line in content: + if line[0] != '#': # skip comments + (x[0],y[0],x[1],y[1]) = map(float,line.split())[8:12] # get start and end coordinates of each segment + boxX[0] = min(boxX[0],x[0],x[1]) + boxX[1] = max(boxX[1],x[0],x[1]) + boxY[0] = min(boxY[0],y[0],y[1]) + boxY[1] = max(boxY[1],y[0],y[1]) + dX = boxX[1]-boxX[0] + dY = boxY[1]-boxY[0] + + scaleImg = imagesize/max(dX,dY) + scalePatch = size/dY + + if scaleImg > 0: # create image + img = Image.new("RGB",map(lambda x:int(round(x))+border*2,(scaleImg*dX,scaleImg*dY)),(255,255,255)) + draw = ImageDraw.Draw(img) + +# read segments and draw them + segment = 0 + connectivityXY = {"0": {"0":[],"%g"%dY:[],},\ + "%g"%dX: {"0":[],"%g"%dY:[],},} + connectivityYX = {"0": {"0":[],"%g"%dX:[],},\ + "%g"%dY: {"0":[],"%g"%dX:[],},} + grainNeighbors = [] + + for line in content: + if line[0] != '#': # skip comments + (x[0],y[0],x[1],y[1]) = map(float,line.split())[8:12] # get start and end coordinates of each segment + # make relative to origin of bounding box + x[0] -= boxX[0] + x[1] -= boxX[0] + y[0]=boxY[1]-y[0] + y[1]=boxY[1]-y[1] + grainNeighbors.append(map(int,line.split()[12:14])) # remember right and left grain per segment + for i in range(2): # store segment to both points + match = False # check whether point is already known (within a small range) + for posX in connectivityXY.keys(): + if (abs(float(posX)-x[i]) 0: + draw.line(map(lambda x:int(scaleImg*x)+border,[x[0],y[0],x[1],y[1]]),fill=(128,128,128)) + draw.text(map(lambda x:int(scaleImg*x)+border,[(x[1]+x[0])/2.0,(y[1]+y[0])/2.0]),"%i"%segment,fill=(0,0,128)) + segment += 1 + + +# top border + keyId = "0" + boundary = connectivityYX[keyId].keys() + boundary.sort(key=float) + for indexBdy in range(len(boundary)-1): + connectivityXY[boundary[indexBdy]][keyId].append(segment) + connectivityXY[boundary[indexBdy+1]][keyId].append(segment) + connectivityYX[keyId][boundary[indexBdy]].append(segment) + connectivityYX[keyId][boundary[indexBdy+1]].append(segment) + if scaleImg > 0: + draw.line(map(lambda x:int(scaleImg*x)+border,[float(boundary[indexBdy]),float(keyId),float(boundary[indexBdy+1]),float(keyId)]),width=3,fill=(128,128*(segment%2),0)) + draw.text(map(lambda x:int(scaleImg*x)+border,[(float(boundary[indexBdy])+float(boundary[indexBdy+1]))/2.0,float(keyId)]),"%i"%segment,fill=(0,0,128)) + segment += 1 + +# right border + keyId = "%g"%(boxX[1]-boxX[0]) + boundary = connectivityXY[keyId].keys() + boundary.sort(key=float) + for indexBdy in range(len(boundary)-1): + connectivityYX[boundary[indexBdy]][keyId].append(segment) + connectivityYX[boundary[indexBdy+1]][keyId].append(segment) + connectivityXY[keyId][boundary[indexBdy]].append(segment) + connectivityXY[keyId][boundary[indexBdy+1]].append(segment) + if scaleImg > 0: + draw.line(map(lambda x:int(scaleImg*x)+border,[float(keyId),float(boundary[indexBdy]),float(keyId),float(boundary[indexBdy+1])]),width=3,fill=(128,128*(segment%2),0)) + draw.text(map(lambda x:int(scaleImg*x)+border,[float(keyId),(float(boundary[indexBdy])+float(boundary[indexBdy+1]))/2.0]),"%i"%segment,fill=(0,0,128)) + segment += 1 + +# bottom border + keyId = "%g"%(boxY[1]-boxY[0]) + boundary = connectivityYX[keyId].keys() + boundary.sort(key=float,reverse=True) + for indexBdy in range(len(boundary)-1): + connectivityXY[boundary[indexBdy]][keyId].append(segment) + connectivityXY[boundary[indexBdy+1]][keyId].append(segment) + connectivityYX[keyId][boundary[indexBdy]].append(segment) + connectivityYX[keyId][boundary[indexBdy+1]].append(segment) + if scaleImg > 0: + draw.line(map(lambda x:int(scaleImg*x)+border,[float(boundary[indexBdy]),float(keyId),float(boundary[indexBdy+1]),float(keyId)]),width=3,fill=(128,128*(segment%2),0)) + draw.text(map(lambda x:int(scaleImg*x)+border,[(float(boundary[indexBdy])+float(boundary[indexBdy+1]))/2.0,float(keyId)]),"%i"%segment,fill=(0,0,128)) + segment += 1 + +# left border + keyId = "0" + boundary = connectivityXY[keyId].keys() + boundary.sort(key=float,reverse=True) + for indexBdy in range(len(boundary)-1): + connectivityYX[boundary[indexBdy]][keyId].append(segment) + connectivityYX[boundary[indexBdy+1]][keyId].append(segment) + connectivityXY[keyId][boundary[indexBdy]].append(segment) + connectivityXY[keyId][boundary[indexBdy+1]].append(segment) + if scaleImg > 0: + draw.line(map(lambda x:int(scaleImg*x)+border,[float(keyId),float(boundary[indexBdy]),float(keyId),float(boundary[indexBdy+1])]),width=3,fill=(128,128*(segment%2),0)) + draw.text(map(lambda x:int(scaleImg*x)+border,[float(keyId),(float(boundary[indexBdy])+float(boundary[indexBdy+1]))/2.0]),"%i"%segment,fill=(0,0,128)) + segment += 1 + + + allkeysX = connectivityXY.keys() + allkeysX.sort() + points = [] + segments = [[] for i in range(segment)] + pointId = 0 + for keyX in allkeysX: + allkeysY = connectivityXY[keyX].keys() + allkeysY.sort() + for keyY in allkeysY: + points.append({'coords': [float(keyX)*scalePatch,float(keyY)*scalePatch], 'segments': connectivityXY[keyX][keyY]}) + for segment in connectivityXY[keyX][keyY]: + if (segments[segment] == None): + segments[segment] = pointId + else: + segments[segment].append(pointId) + if scaleImg > 0: + draw.text(map(lambda x:int(scaleImg*x)+border,[float(keyX),float(keyY)]),"%i"%pointId,fill=(0,0,0)) + pointId += 1 + + if scaleImg > 0: + img.save(imagename+'.png',"PNG") + + + grains = {'draw': [], 'legs': []} + pointId = 0 + for point in points: + while point['segments']: + myStart = pointId + grainDraw = [points[myStart]['coords']] + innerAngleSum = 0.0 + myWalk = point['segments'].pop() + grainLegs = [myWalk] + if segments[myWalk][0] == myStart: + myEnd = segments[myWalk][1] + else: + myEnd = segments[myWalk][0] + + while (myEnd != pointId): + myV = [points[myEnd]['coords'][0]-points[myStart]['coords'][0],\ + points[myEnd]['coords'][1]-points[myStart]['coords'][1]] + myLen = math.sqrt(myV[0]**2+myV[1]**2) + best = {'product': -2.0, 'peek': -1, 'len': -1, 'point': -1} + for peek in points[myEnd]['segments']: + if peek == myWalk: + continue + if segments[peek][0] == myEnd: + peekEnd = segments[peek][1] + else: + peekEnd = segments[peek][0] + peekV = [points[myEnd]['coords'][0]-points[peekEnd]['coords'][0],\ + points[myEnd]['coords'][1]-points[peekEnd]['coords'][1]] + peekLen = math.sqrt(peekV[0]**2+peekV[1]**2) + crossproduct = (myV[0]*peekV[1]-myV[1]*peekV[0])/myLen/peekLen + dotproduct = (myV[0]*peekV[0]+myV[1]*peekV[1])/myLen/peekLen + if crossproduct*(dotproduct+1.0) >= best['product']: + best['product'] = crossproduct*(dotproduct+1.0) + best['peek'] = peek + best['point'] = peekEnd + innerAngleSum += best['product'] + myWalk = best['peek'] + myStart = myEnd + myEnd = best['point'] + points[myStart]['segments'].remove(myWalk) + grainDraw.append(points[myStart]['coords']) + grainLegs.append(myWalk) + if innerAngleSum > 0.0: + grains['draw'].append(grainDraw) + grains['legs'].append(grainLegs) + else: + grains['box'] = grainLegs + pointId += 1 + + +# build overall data structure + + rcData = {'dimension':[dX,dY], 'point': [],'segment': [], 'grain': [], 'grainMapping': []} + + for point in points: + rcData['point'].append(point['coords']) + print "found %i points"%(len(rcData['point'])) + + for segment in segments: + rcData['segment'].append(segment) + print "built %i segments"%(len(rcData['segment'])) + + for grain in grains['legs']: + rcData['grain'].append(grain) + myNeighbors = {} + for leg in grain: + if leg < len(grainNeighbors): + for side in range(2): + if grainNeighbors[leg][side] in myNeighbors: + myNeighbors[grainNeighbors[leg][side]] += 1 + else: + myNeighbors[grainNeighbors[leg][side]] = 1 + if myNeighbors: # do I have any neighbors + rcData['grainMapping'].append(sorted(myNeighbors.iteritems(), key=lambda (k,v): (v,k), reverse=True)[0][0]) # most frequent grain is me + print "found %i grains"%(len(rcData['grain'])) + + rcData['box'] = grains['box'] + + + if scaleImg > 0: + grainID = 0 + for grain in grains['draw']: + coords = [0,0] + for point in grain: + coords[0] += int(scaleImg/scalePatch*point[0]) + coords[1] += int(scaleImg/scalePatch*point[1]) + coords[0] /= len(grain) + coords[1] /= len(grain) + draw.text(map(lambda x:x+border,coords),'%i -> %i'%(grainID,rcData['grainMapping'][grainID]),fill=(128,32,32)) + grainID += 1 + + img.save(os.path.splitext(args[0])[0]+'.png',"PNG") + + return rcData + + +def init(): + return ["*new_model yes", + "*select_clear", + "*reset", + "*set_nodes off", + "*elements_solid", + "*show_view 4", + "*reset_view", + "*view_perspective", + "*redraw", + ] + + +def sample(a,n,size): + + cmds = [\ +# gauge + "*add_points %f %f %f"%(-size*a,size*a,0), + "*add_points %f %f %f"%( size*a,size*a,0), + "*add_points %f %f %f"%( size*a,-size*a,0), + "*add_points %f %f %f"%(-size*a,-size*a,0), + "*set_curve_type line", + "*add_curves %i %i"%(1,2), + "*add_curves %i %i"%(3,4), + "*set_curve_div_type_fix_ndiv", + "*set_curve_div_num %i"%n, + "*apply_curve_divisions", + "1 2 #", + "*add_curves %i %i"%(2,3), # right side + "*add_curves %i %i"%(4,1), # left side + "*set_curve_div_type_fix_ndiv", + "*set_curve_div_num %i"%n, + "*apply_curve_divisions", + "3 4 #", + ] + + return cmds + + +def patch(a,n,mesh,rcData): + cmds = [] + for l in range(len(rcData['point'])): # generate all points + cmds.append("*add_points %f %f %f"%(rcData['point'][l][0]-a*rcData['dimension'][0]/rcData['dimension'][1]/2.0,rcData['point'][l][1]-a/2.0,0)) + + cmds.append(["*set_curve_type line", + "*set_curve_div_type_fix_ndiv", + ]) + for m in range(len(rcData['segment'])): # generate all curves and subdivide them for overall balanced piece length + start = rcData['segment'][m][0] + end = rcData['segment'][m][1] + cmds.append([\ + "*add_curves %i %i" %(start+rcData['offsetPoints'], + end +rcData['offsetPoints']), + "*set_curve_div_num %i"%(max(1,round(math.sqrt((rcData['point'][start][0]-rcData['point'][end][0])**2+\ + (rcData['point'][start][1]-rcData['point'][end][1])**2)/a*n))), + "*apply_curve_divisions", + "%i #"%(m+rcData['offsetSegments']), + ]) + + grain = 0 + cmds.append('(!)locals["last"] = py_get_int("nelements()")') + for g in rcData['grain']: + cmds.append([\ + '(!)locals["first"] = locals["last"]+1', + "*%s "%mesh+" ".join([str(rcData['offsetSegments']+x) for x in g])+" #", + '(!)locals["last"] = py_get_int("nelements()")', + "*select_elements", + '(?)"%i to %i #"%(locals["first"],locals["last"])', + "*store_elements grain_%i"%rcData['grainMapping'][grain], + "all_selected", + "*select_clear", + ]) + grain += 1 + + return cmds + + +def gage(mesh,rcData): + + return([\ + "*%s "%mesh + + " ".join([str(x) for x in range(1,rcData['offsetSegments'])]) + + " " + + " ".join([str(rcData['offsetSegments']+x)for x in rcData['box']]) + + " #", + "*select_reset", + "*select_clear", + "*select_elements", + "all_existing", + "*select_mode_except", + ['grain_%i'%(i+1) for i in range(len(rcData['grain']))], + "#", + "*store_elements matrix", + "all_selected", + "*select_mode_invert", + "*select_elements", + "all_existing", + "*store_elements grains", + "all_selected", + "*select_clear", + "*select_reset", + ]) + + +def expand3D(thickness,steps): + return([\ + "*set_expand_translation z %f"%(thickness/steps), + "*set_expand_repetitions %i"%steps, + "*expand_elements", + "all_existing", + ]) + + +def initial_conditions(grainNumber,grainMapping): + cmds = [\ + "*new_icond", + "*icond_name temperature", + "*icond_type state_variable", + "*icond_param_value state_var_id 1", + "*icond_dof_value var 300", + "*add_icond_elements", + "all_existing", + "*new_icond", + "*icond_name homogenization", + "*icond_type state_variable", + "*icond_param_value state_var_id 2", + "*icond_dof_value var 1", + "*add_icond_elements", + "all_existing", + ] + + for grain in range(grainNumber): + cmds.append([\ + "*new_icond", + "*icond_name grain_%i"%grainMapping[grain], + "*icond_type state_variable", + "*icond_param_value state_var_id 3", + "*icond_dof_value var %i"%(grain+1), + "*add_icond_elements", + "grain_%i"%grainMapping[grain], + "", + ]) + cmds.append([\ + "*new_icond", + "*icond_name rim", + "*icond_type state_variable", + "*icond_param_value state_var_id 3", + "*icond_dof_value var %i"%(grainNumber+1), + "*add_icond_elements", + "matrix", + ]) + return cmds + + +def boundary_conditions(rate,thickness, a,size): + inner = size*(1 - 1.0e-4) * a + outer = size*(1 + 1.0e-4) * a + upper = size*(1 + 1.0e-4) * a + lower = size*(1 - 1.0e-4) * a + + return [\ + "*new_md_table 1 1", + "*table_name linear", + "*set_md_table_type 1 time", + "*table_add", + "0 0", + "1 1", + "*select_method_box", + "*new_apply", + "*apply_name pull_bottom", + "*apply_type fixed_displacement", + "*apply_dof y", + "*apply_dof_value y %f"%(-rate*(inner+outer)/2.0), + "*apply_dof_table y linear", + "*select_clear_nodes", + "*select_nodes", + "%f %f"%(-outer,outer), + "%f %f"%(-upper,-lower), + "%f %f"%(-.0001*a,(thickness+1.0e-4)*a), + "*add_apply_nodes", + "all_selected", + "*new_apply", + "*apply_name pull_top", + "*apply_type fixed_displacement", + "*apply_dof y", + "*apply_dof_value y %f"%(rate*(inner+outer)/2.0), + "*apply_dof_table y linear", + "*select_clear_nodes", + "*select_nodes", + "%f %f"%(-outer,outer), + "%f %f"%(lower,upper), + "%f %f"%(-.0001*a,(thickness+1.0e-4)*a), + "*add_apply_nodes", + "all_selected", + "*new_apply", + "*apply_name fix_x", + "*apply_type fixed_displacement", + "*apply_dof x", + "*apply_dof_value x 0", + "*select_clear_nodes", + "*select_nodes", + "%f %f"%(-outer,-inner), + "%f %f"%(lower,upper), + "%f %f"%(-.0001*a,.0001*a), + "%f %f"%(-outer,-inner), + "%f %f"%(lower,upper), + "%f %f"%((thickness-1.0e-4)*a,(thickness+1.0e-4)*a), + "%f %f"%(-outer,-inner), + "%f %f"%(-upper,-lower), + "%f %f"%(-.0001*a,.0001*a), + "%f %f"%(-outer,-inner), + "%f %f"%(-upper,-lower), + "%f %f"%((thickness-1.0e-4)*a,(thickness+1.0e-4)*a), + "*add_apply_nodes", + "all_selected", + "*new_apply", + "*apply_name fix_z", + "*apply_type fixed_displacement", + "*apply_dof z", + "*apply_dof_value z 0", + "*select_clear_nodes", + "*select_nodes", + "%f %f"%(-outer,-inner), + "%f %f"%(lower,upper), + "%f %f"%(-.0001*a,.0001*a), + "%f %f"%(-outer,-inner), + "%f %f"%(-upper,-lower), + "%f %f"%(-.0001*a,.0001*a), + "%f %f"%(inner,outer), + "%f %f"%(lower,upper), + "%f %f"%(-.0001*a,.0001*a), + "%f %f"%(inner,outer), + "%f %f"%(-upper,-lower), + "%f %f"%(-.0001*a,.0001*a), + "*add_apply_nodes", + "all_selected", + "*select_clear", + "*select_reset", + ] + +def materials(): + return [\ + "*new_material", + "*material_name patch", + "*material_type mechanical:hypoelastic", + "*material_option hypoelastic:method:hypela2", + "*material_option hypoelastic:pass:def_rot", + "*add_material_elements", + "all_existing", + ] + + +def loadcase(time,incs,Ftol): + return [\ + "*new_loadcase", + "*loadcase_name puller", + "*loadcase_type static", + "*loadcase_value time", + "%g"%time, + "*loadcase_value nsteps", + "%i"%incs, + "*loadcase_value maxrec", + "20", + "*loadcase_value ntime_cuts", + "30", + "*loadcase_value force", + "%g"%Ftol, + ] + + +def job(grainNumber,grainMapping): + return [\ + "*new_job", + "*job_name pull", + "*job_class mechanical", + "*add_job_loadcases puller", + "*add_job_iconds homogenization", + ["*add_job_iconds grain_%i"%i for i in grainMapping[:grainNumber]], + "*add_job_iconds rim", + "*job_option dimen:three | 3D analysis", + "*job_option strain:large | finite strains", + "*job_option large_strn_proc:upd_lagrange | updated Lagrange framework", + "*job_option plas_proc:multiplicative | multiplicative decomp of F", + "*job_option solver_nonsym:on | nonsymmetrical solution", + "*job_option solver:mfront_sparse | multi-frontal sparse", + "*job_param stef_boltz 5.670400e-8", + "*job_param univ_gas_const 8.314472", + "*job_param planck_radiation_2 1.4387752e-2", + "*job_param speed_light_vacuum 299792458", + "*job_usersub_file /san/%s/FEM/subroutine_svn/mpie_cpfem_marc2007r1.f90 | subroutine definition"%(pwd.getpwuid(os.geteuid())[0].rpartition("\\")[2]), + "*job_option user_source:compile_save", + ] + + +def postprocess(): + return [\ + "*add_post_tensor stress", + "*add_post_tensor strain", + "*add_post_var von_mises", + "", + ] + + + +def cleanUp(a): + return [\ + "*remove_curves", + "all_existing", + "*remove_points", + "all_existing", + "*set_sweep_tolerance %f"%(1e-3*a), + "*sweep_all", + "*renumber_all", + ] + + +# ----------------------- MAIN ------------------------------- + +parser = OptionParser() +parser.add_option("-p", "--port", type="int",\ + dest="port",\ + help="Mentat connection port") +parser.add_option("-a", "--patchsize", type="float",\ + dest="size",\ + help="height of patch [%default]") +parser.add_option("-f", "--frame", type="float",\ + dest="frame",\ + help="frame thickness in units of patch height [%default]") +parser.add_option("-n", "--resolution", type="int",\ + dest="resolution",\ + help="number of elements along patch perimeter [%default]") +parser.add_option("-e", "--strain", type="float",\ + dest="strain",\ + help="final strain to reach in simulation [%default]") +parser.add_option("-r", "--rate", type="float",\ + dest="strainrate",\ + help="(engineering) strain rate to simulate") +parser.add_option("-i", "--increments", type="int",\ + dest="increments",\ + help="number of increments to take") +parser.add_option("-s", "--imagesize", type="int",\ + dest="imgsize",\ + help="size of image") +parser.add_option("-b", "--border", type="int",\ + dest="border",\ + help="border of image") +parser.add_option("-t", "--tolerance", type="float",\ + dest="tolerance",\ + help="relative tolerance of pixel positions to be swept") +parser.add_option("-m", "--mesh", choices=['dt_planar_trimesh','af_planar_trimesh','af_planar_quadmesh'],\ + dest="mesh",\ + help="algorithm and element type for automeshing [%default]") + +parser.set_defaults(size = 1.0) +parser.set_defaults(frame = 0.5) +parser.set_defaults(resolution = 30) +parser.set_defaults(strain = 0.2) +parser.set_defaults(strainrate = 1.0e-3) +parser.set_defaults(increments = 200) +parser.set_defaults(imgsize = 0) +parser.set_defaults(border = 20) +parser.set_defaults(tolerance = 1.0e-3) +parser.set_defaults(mesh = 'dt_planar_trimesh') + +(options, args) = parser.parse_args() +if not len(args): + parser.error('no boundary file specified') + +try: + boundaryFile = open(args[0]) + boundarySegments = boundaryFile.readlines() + boundaryFile.close() +except: + print 'unable to read boundary file "%s"'%args[0] + sys.exit(-1) + +myName = os.path.splitext(args[0])[0] +print "\n",myName +orientationData = rcbOrientationParser(boundarySegments) +rcData = rcbParser(boundarySegments,options.size,options.tolerance,myName,options.imgsize,options.border) + +# ----- write texture data to file ----- +configFile = open(os.path.splitext(args[0])[0]+'.config','w') +configFile.write('\n\n\n\n') +for i,grain in enumerate(rcData['grainMapping']): + configFile.write('\n[grain %i]\n'%grain) + configFile.write('(constituent)\tphase %i\ttexture %i\tfraction 1.0\n'%(i+1,i+1)) + +configFile.write('\n\n\n\n') +for grain in rcData['grainMapping']: + configFile.write('\n[grain %i]\n'%grain) + +configFile.write('\n\n\n\n') +for grain in rcData['grainMapping']: + configFile.write('\n[grain %i]\n'%grain) + configFile.write('(gauss)\tphi1\t%f\tphi\t%f\tphi2\t%f\tscatter\t0.0\tfraction\t1.0\n'\ + %(math.degrees(orientationData[grain-1][0]),math.degrees(orientationData[grain-1][1]),math.degrees(orientationData[grain-1][2]))) +configFile.close() + +rcData['offsetPoints'] = 1+4 # gage definition generates 4 points +rcData['offsetSegments'] = 1+4 # gage definition generates 4 segments + +cmds = [\ + init(), + sample(options.size,12,options.frame+0.5), + patch(options.size,options.resolution,options.mesh,rcData), + gage(options.mesh,rcData), + expand3D(options.size/6,4), + cleanUp(options.size), + materials(), + initial_conditions(len(rcData['grain']),rcData['grainMapping']), + boundary_conditions(options.strainrate,options.size/6,options.size,options.frame+0.5), + loadcase(options.strain/options.strainrate,options.increments,0.03), + job(len(rcData['grain']),rcData['grainMapping']), + postprocess(), + ["*identify_sets","*regen","*fill_view","*save_as_model %s yes"%(myName)], +] + +outputLocals = {} +if (options.port != None): + py_connect('',options.port) + output(cmds,outputLocals,'Mentat') + py_disconnect() +else: + output(cmds,outputLocals,'Stdout') + +print outputLocals + +# "*job_option large:on | large displacement", +# "*job_option plasticity:l_strn_mn_add | large strain additive", +# "*job_option cdilatation:on | constant dilatation", +# "*job_option update:on | updated lagrange procedure", +# "*job_option finite:on | large strains", +# "*job_option restart_mode:write | enable restarting", diff --git a/processing/pre/mentat_pbcOnBoxMesh b/processing/pre/mentat_pbcOnBoxMesh new file mode 100755 index 000000000..5fda5adfb --- /dev/null +++ b/processing/pre/mentat_pbcOnBoxMesh @@ -0,0 +1,175 @@ +#!/usr/bin/env python + +releases = ['2010','2008r1','2007r1','2005r3'] + + +import sys,os,pwd,math,re +#import Image,ImageDraw +from optparse import OptionParser + +for release in releases: + libPath = '/msc/mentat%s/shlib/'%release + if os.path.exists(libPath): + sys.path.append(libPath) + for subdir in [os.path.join(libPath,file) + for file in os.listdir(libPath) + if os.path.isdir(os.path.join(libPath,file))]: + sys.path.append(subdir) + break +from py_mentat import * + +def outMentat(cmd,locals): + if cmd[0:3] == '(!)': + exec(cmd[3:]) + elif cmd[0:3] == '(?)': + cmd = eval(cmd[3:]) + py_send(cmd) + else: + py_send(cmd) + return + +def outStdout(cmd,locals): + if cmd[0:3] == '(!)': + exec(cmd[3:]) + elif cmd[0:3] == '(?)': + cmd = eval(cmd[3:]) + print cmd + else: + print cmd + return + + +def output(cmds,locals,dest): + for cmd in cmds: + if isinstance(cmd,list): + output(cmd,locals,dest) + else: + {\ + 'Mentat': outMentat,\ + 'Stdout': outStdout,\ + }[dest](cmd,locals) + return + + + +def servoLink(): + + cmds = [] + base = ['x','y','z'] + box = {'min': {'x': float(sys.maxint),'y': float(sys.maxint),'z': float(sys.maxint)}, + 'max': {'x':-float(sys.maxint),'y':-float(sys.maxint),'z':-float(sys.maxint)}, + 'delta': {'x':0,'y':0,'z':0}, + } + Nnodes = py_get_int("nnodes()") + NodeCoords = [{'x':py_get_float("node_x(%i)"%(node)), + 'y':py_get_float("node_y(%i)"%(node)), + 'z':py_get_float("node_z(%i)"%(node)),} for node in range(1,1+Nnodes)] + + for node in range(Nnodes): # find the bounding box + for coord in base: # check each direction in turn + box['min'][coord] = min(box['min'][coord],NodeCoords[node][coord]) + box['max'][coord] = max(box['max'][coord],NodeCoords[node][coord]) + + for coord in base: # calc the dimension of the bounding box + box['delta'][coord] = box['max'][coord] - box['min'][coord] + + baseNode = {} + linkNodes = [] + + for node in range(Nnodes): # loop over all nodes + pos = {} + key = {} + maxFlag = {'x': False, 'y': False, 'z': False} + Nmax = 0 + Nmin = 0 + for coord in base: # for each direction + key[coord] = "%.8e"%NodeCoords[node][coord] # translate position to string + if (key[coord] == "%.8e"%box['min'][coord]): # compare to min of bounding box (i.e. is on outer face?) + Nmin += 1 # count outer (back) face membership + elif (key[coord] == "%.8e"%box['max'][coord]): # compare to max of bounding box (i.e. is on outer face?) + Nmax += 1 # count outer (front) face memebership + maxFlag[coord] = True # remember face membership (for linked nodes) + + if Nmin > 0 and Nmin > Nmax: # node is on more back than font faces + # prepare for any non-existing entries in the data structure + if key['x'] not in baseNode.keys(): + baseNode[key['x']] = {} + if key['y'] not in baseNode[key['x']].keys(): + baseNode[key['x']][key['y']] = {} + if key['z'] not in baseNode[key['x']][key['y']].keys(): + baseNode[key['x']][key['y']][key['z']] = 0 + + baseNode[key['x']][key['y']][key['z']] = node+1 # remember the base node id + + elif Nmax > 0 and Nmax >= Nmin: # node is on at least as many front than back faces + linkNodes.append({'id': node+1,'coord': NodeCoords[node], 'onFaces': Nmax,'faceMember': maxFlag}) + + + baseCorner = baseNode["%.8e"%box['min']['x']]["%.8e"%box['min']['y']]["%.8e"%box['min']['z']] # detect ultimate base node + + for node in linkNodes: # loop over all linked nodes + linkCoord = [node['coord']] # start list of control node coords with my coords + for dir in base: # check for each direction + if node['faceMember'][dir]: # me on this front face + linkCoord[0][dir] = box['min'][dir] # project me onto rear face along dir + linkCoord.append({'x':box['min']['x'],'y':box['min']['y'],'z':box['min']['z'],}) # append base corner + linkCoord[-1][dir] = box['max'][dir] # stretch it to corresponding control leg of "dir" + + nLinks = len(linkCoord) + + for dof in [1,2,3]: + cmds.append([ + "*new_link *link_class servo", + "*link_class servo *tied_node %i"%node['id'], + "*link_class servo *tied_dof %i"%dof, + "*servo_nterms %i"%(1+nLinks), + ]) + for i in range(nLinks): + cmds.append([ + "*link_class servo *servo_ret_node %i %i"%(i+1,baseNode["%.8e"%linkCoord[i]['x']]["%.8e"%linkCoord[i]['y']]["%.8e"%linkCoord[i]['z']]), + "*link_class servo *servo_ret_dof %i %i"%(i+1,dof), + "*link_class servo *servo_ret_coef %i 1"%(i+1), + ]) + cmds.append([ + "*link_class servo *servo_ret_node %i %i"%(1+nLinks,baseCorner), + "*link_class servo *servo_ret_dof %i %i"%(1+nLinks,dof), + "*link_class servo *servo_ret_coef %i -%i"%(1+nLinks,nLinks-1), + ]) + + cmds.append([ + "*select_nodes", + ["%i"%node['id'] for node in linkNodes], + "#", + ]) + + return cmds + + + +# ----------------------- MAIN ------------------------------- + +parser = OptionParser() +parser.add_option("-p", "--port", type="int",\ + dest="port",\ + help="Mentat connection port [%default]") +parser.add_option("-v", "--verbose", action="store_true",\ + dest="verbose",\ + help="write Mentat command stream also to stdout [%default]") +parser.set_defaults(port=40007) +parser.set_defaults(verbose=False) + +(options, args) = parser.parse_args() + +outputLocals = {} +print 'waiting to connect...' +py_connect('',options.port) +output(['*remove_all_servos', + '*sweep_all', + '*renumber_nodes'],outputLocals,'Mentat') # script depends on consecutive numbering of nodes +cmds = servoLink() +print 'connected...' +output(cmds,outputLocals,'Mentat') +py_disconnect() + +if options.verbose: + output(cmds,outputLocals,'Stdout')