DAMASK_EICMD/processing/pre/patchFromReconstructedBound...

982 lines
37 KiB
Python
Executable File

#!/usr/bin/env python
import sys,os,math,re,string, damask
from optparse import OptionParser, OptionGroup, Option, SUPPRESS_HELP
try: # check for Python Image Lib
import Image,ImageDraw
ImageCapability = True
except:
ImageCapability = False
sys.path.append(damask.solver.Marc().libraryPath('../../'))
try: # check for MSC.Mentat Python interface
from py_mentat import *
MentatCapability = True
except:
MentatCapability = False
# -----------------------------
class extendedOption(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)
def outMentat(cmd,locals):
if cmd[0:3] == '(!)':
exec(cmd[3:])
elif cmd[0:3] == '(?)':
cmd = eval(cmd[3:])
py_send(cmd)
if 'log' in locals: locals['log'].append(cmd)
else:
py_send(cmd)
if 'log' in locals: locals['log'].append(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:
m = re.match(r'\s*(#|$)',line)
if m: continue # skip comments and blank lines
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,M,size,tolerance): # 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:
m = re.match(r'\s*(#|$)',line)
if m: continue # skip comments and blank lines
(x[0],y[0],x[1],y[1]) = map(float,line.split())[8:12] # get start and end coordinates of each segment.
(x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0]) # apply transformation to coordinates
(x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1]) # to get rcb --> Euler system
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]
scalePatch = size/dX
# read segments
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:
m = re.match(r'\s*(#|$)',line)
if m: continue # skip comments and blank lines
(x[0],y[0],x[1],y[1]) = map(float,line.split())[8:12] # get start and end coordinates of each segment
(x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0]) # apply transformation to coordinates
(x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1]) # to get rcb --> Euler system
x[0] -= boxX[0] # make relative to origin of bounding box
x[1] -= boxX[0]
y[0] -= boxY[0]
y[1] -= boxY[0]
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])<dX*tolerance):
for posY in connectivityXY[posX].keys():
if (abs(float(posY)-y[i])<dY*tolerance):
keyX = posX
keyY = posY
match = True
break
break
if (not match):
# force to boundary if inside tolerance to it
if (abs(x[i])<dX*tolerance):
x[i] = 0
if (abs(dX-x[i])<dX*tolerance):
x[i] = dX
if (abs(y[i])<dY*tolerance):
y[i] = 0
if (abs(dY-y[i])<dY*tolerance):
y[i] = dY
keyX = "%g"%x[i]
keyY = "%g"%y[i]
if keyX not in connectivityXY: # create new hash entry for so far unknown point
connectivityXY[keyX] = {}
if keyY not in connectivityXY[keyX]: # create new hash entry for so far unknown point
connectivityXY[keyX][keyY] = []
if keyY not in connectivityYX: # create new hash entry for so far unknown point
connectivityYX[keyY] = {}
if keyX not in connectivityYX[keyY]: # create new hash entry for so far unknown point
connectivityYX[keyY][keyX] = []
connectivityXY[keyX][keyY].append(segment)
connectivityYX[keyY][keyX].append(segment)
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)
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)
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)
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)
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)
pointId += 1
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']
if myWalk in points[myStart]['segments']:
points[myStart]['segments'].remove(myWalk)
else:
sys.stderr.write(str(myWalk)+' not in segments of '+str(myStart))
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 legs in grains['legs']: # loop over grains
rcData['grain'].append(legs) # store list of boundary segments
myNeighbors = {}
for leg in legs: # test each boundary segment
if leg < len(grainNeighbors): # a valid segment index?
for side in range(2): # look at both sides of the segment
if grainNeighbors[leg][side] in myNeighbors: # count occurrence of grain IDs
myNeighbors[grainNeighbors[leg][side]] += 1
else:
myNeighbors[grainNeighbors[leg][side]] = 1
if myNeighbors: # do I have any neighbors (i.e., non-bounding box segment)
candidateGrains = sorted(myNeighbors.iteritems(), key=lambda (k,v): (v,k), reverse=True) # sort grain counting
if candidateGrains[0][0] not in rcData['grainMapping']: # most frequent one not yet seen?
rcData['grainMapping'].append(candidateGrains[0][0]) # must be me then
else:
rcData['grainMapping'].append(candidateGrains[1][0]) # special case of bi-crystal situation...
print " found %i grains\n"%(len(rcData['grain']))
rcData['box'] = grains['box']
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(size,aspect,n,xmargin,ymargin):
cmds = [\
# gauge
"*add_points %f %f %f"%(-size*(0.5+xmargin), size*(0.5*aspect+ymargin),0),
"*add_points %f %f %f"%( size*(0.5+xmargin), size*(0.5*aspect+ymargin),0),
"*add_points %f %f %f"%( size*(0.5+xmargin),-size*(0.5*aspect+ymargin),0),
"*add_points %f %f %f"%(-size*(0.5+xmargin),-size*(0.5*aspect+ymargin),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/2.0,rcData['point'][l][1]-a/rcData['dimension'][0]*rcData['dimension'][1]/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'%rcData['grainMapping'][i] 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, size,aspect,xmargin,ymargin):
inner = (1 - 1.0e-4) * size*(0.5+xmargin)
outer = (1 + 1.0e-4) * size*(0.5+xmargin)
lower = (1 - 1.0e-4) * size*(0.5*aspect+ymargin)
upper = (1 + 1.0e-4) * size*(0.5*aspect+ymargin)
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*(lower+upper)/2.0),
"*apply_dof_table y linear",
"*select_clear_nodes",
"*select_nodes",
"%f %f"%(-outer,outer),
"%f %f"%(-upper,-lower),
"%f %f"%(-.0001*thickness,1.0001*thickness),
"*add_apply_nodes",
"all_selected",
"*new_apply",
"*apply_name pull_top",
"*apply_type fixed_displacement",
"*apply_dof y",
"*apply_dof_value y %f"%(rate*(lower+upper)/2.0),
"*apply_dof_table y linear",
"*select_clear_nodes",
"*select_nodes",
"%f %f"%(-outer,outer),
"%f %f"%(lower,upper),
"%f %f"%(-.0001*thickness,1.0001*thickness),
"*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*thickness,.0001*thickness),
"%f %f"%(-outer,-inner),
"%f %f"%(lower,upper),
"%f %f"%(0.9999*thickness,1.0001*thickness),
"%f %f"%(-outer,-inner),
"%f %f"%(-upper,-lower),
"%f %f"%(-.0001*thickness,.0001*thickness),
"%f %f"%(-outer,-inner),
"%f %f"%(-upper,-lower),
"%f %f"%(0.9999*thickness,1.0001*thickness),
"*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*thickness,.0001*thickness),
"%f %f"%(-outer,-inner),
"%f %f"%(-upper,-lower),
"%f %f"%(-.0001*thickness,.0001*thickness),
"%f %f"%(inner,outer),
"%f %f"%(lower,upper),
"%f %f"%(-.0001*thickness,.0001*thickness),
"%f %f"%(inner,outer),
"%f %f"%(-upper,-lower),
"%f %f"%(-.0001*thickness,.0001*thickness),
"*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,twoD):
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:%s | analysis dimension"%({True:'two ',False:'three'}[twoD]),
"*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/DAMASK/code/mpie_cpfem_marc2010.f90 | subroutine definition"%(pwd.getpwuid(os.geteuid())[0].rpartition("\\")[2]),
"*job_option user_source:compile_save",
]
# "*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",
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-5*a),
"*sweep_all",
"*renumber_all",
]
# -------------------------
def image(name,imgsize,marginX,marginY,rcData):
# -------------------------
dX = max([coords[0] for coords in rcData['point']])
dY = max([coords[1] for coords in rcData['point']])
offsetX = imgsize*marginX
offsetY = imgsize*marginY
sizeX = int(imgsize*(1 +2*marginX))
sizeY = int(imgsize*(dY/dX+2*marginY))
scaleImg = imgsize/dX # rescale from max x coord
img = Image.new("RGB",(sizeX,sizeY),(255,255,255))
draw = ImageDraw.Draw(img)
for id,point in enumerate(rcData['point']):
draw.text([offsetX+point[0]*scaleImg,sizeY-(offsetY+point[1]*scaleImg)],"%i"%id,fill=(0,0,0))
for id,vertex in enumerate(rcData['segment']):
start = rcData['point'][vertex[0]]
end = rcData['point'][vertex[1]]
draw.text([offsetX+(start[0]+end[0])/2.0*scaleImg,sizeY-(offsetY+(start[1]+end[1])/2.0*scaleImg)],"%i"%id,fill=(0,0,128))
draw.line([offsetX+start[0]*scaleImg,sizeY-(offsetY+start[1]*scaleImg),
offsetX+ end[0]*scaleImg,sizeY-(offsetY+ end[1]*scaleImg)],width=1,fill=(128,128,128))
for id,segment in enumerate(rcData['box']):
start = rcData['point'][rcData['segment'][segment][0]]
end = rcData['point'][rcData['segment'][segment][1]]
draw.line([offsetX+start[0]*scaleImg,sizeY-(offsetY+start[1]*scaleImg),
offsetX+ end[0]*scaleImg,sizeY-(offsetY+ end[1]*scaleImg)],width=3,fill=(128,128*(id%2),0))
for grain,origGrain in enumerate(rcData['grainMapping']):
center = [0.0,0.0]
for segment in rcData['grain'][grain]: # loop thru segments around grain
for point in rcData['segment'][segment]: # take start and end points
center[0] += rcData['point'][point][0] # build vector sum
center[1] += rcData['point'][point][1]
center[0] /= len(rcData['grain'][grain])*2 # normalize by two times segment count, i.e. point count
center[1] /= len(rcData['grain'][grain])*2
draw.text([offsetX+center[0]*scaleImg,sizeY-(offsetY+center[1]*scaleImg)],'%i -> %i'%(grain,origGrain),fill=(128,32,32))
img.save(name+'.png',"PNG")
# -------------------------
def inside(x,y,points): # tests whether point(x,y) is within polygon described by points
# -------------------------
inside = False
npoints=len(points)
(x1,y1) = points[npoints-1] # start with last point of points
startover = (y1 >= y) # am I above testpoint?
for i in range(npoints): # loop through all points
(x2,y2) = points[i] # next point
endover = (y2 >= y) # am I above testpoint?
if (startover != endover): # one above one below testpoint?
if((y2 - y)*(x2 - x1) <= (y2 - y1)*(x2 - x)): # check for intersection
if (endover):
inside = not inside # found intersection
else:
if (not endover):
inside = not inside # found intersection
startover = endover # make second point first point
(x1,y1) = (x2,y2)
return inside
# -------------------------
def fftbuild(rcData, height,xframe,yframe,resolution,extrusion): # build array of grain numbers
# -------------------------
maxX = -1.*sys.maxint
maxY = -1.*sys.maxint
for line in rcData['point']: # find data range
(x,y) = line
maxX = max(maxX, x)
maxY = max(maxY, y)
xsize = maxX+2*xframe # add framsize
ysize = maxY+2*yframe
xres=round(resolution/2.0)*2 # use only even resolution
yres=round(xres/xsize*ysize/2.0)*2 # calculate other resolutions
zres = 1
zsize = min([xsize/xres,ysize/yres])
fftdata = {'fftpoints':[], \
'resolution':(xres,yres,zres), \
'dimension':(xsize,ysize,zsize)}
frameindex=len(rcData['grain'])+1 # calculate frame index as largest grain index plus one
dx = xsize/(xres+1) # calculate step sizes
dy = ysize/(yres+1)
grainpoints = []
for segments in rcData['grain']: # get segments of each grain
points = {}
for i,segment in enumerate(segments[:-1]): # loop thru segments except last (s=[start,end])
points[rcData['segment'][segment][0]] = i # assign segment index to start point
points[rcData['segment'][segment][1]] = i # assigne segment index to endpoint
for i in range(2): # check points of last segment
if points[rcData['segment'][segments[-1]][i]] != 0: # not on first segment
points[rcData['segment'][segments[-1]][i]] = len(segments)-1 # assign segment index to last point
grainpoints.append([]) # start out blank for current grain
for p in sorted(points, key=points.get): # loop thru set of sorted points
grainpoints[-1].append([rcData['point'][p][0],rcData['point'][p][1]]) # append x,y of point
bestGuess = 0 # assume grain 0 as best guess
for i in range(int(xres*yres)): # walk through all points in xy plane
xtest = -xframe+((i%xres)+0.5)*dx # calculate coordinates
ytest = -yframe+(int(i/xres)+0.5)*dy
if(xtest < 0 or xtest > maxX): # check wether part of frame
if( ytest < 0 or ytest > maxY): # part of edges
fftdata['fftpoints'].append(frameindex+2) # append frameindex to result array
else: # part of xframe
fftdata['fftpoints'].append(frameindex) # append frameindex to result array
elif( ytest < 0 or ytest > maxY): # part of yframe
fftdata['fftpoints'].append(frameindex+1) # append frameindex to result array
else:
if inside(xtest,ytest,grainpoints[bestGuess]): # check best guess first
fftdata['fftpoints'].append(bestGuess+1)
else: # no success
for g in range(len(grainpoints)): # test all
if inside(xtest,ytest,grainpoints[g]):
fftdata['fftpoints'].append(g+1)
bestGuess = g
break
return fftdata
# ----------------------- MAIN -------------------------------
parser = OptionParser(option_class=extendedOption, usage='%prog [options] datafile[s]', description = """
Produce image, spectral geometry description, and (auto) Mentat procedure from TSL/OIM
reconstructed boundary file
""" + string.replace('$Id$','\n','\\n')
)
parser.add_option("-o", "--output", action='extend', dest='output', type='string', \
help="types of output [image,mentat,procedure,spectral]")
parser.add_option("-p", "--port", type="int",\
dest="port",\
help="Mentat connection port")
parser.add_option("-2", "--twodimensional", action="store_true",\
dest="twoD",\
help="twodimensional model [%default]")
parser.add_option("-s","--patchsize", type="float",\
dest="size",\
help="height of patch [%default]")
parser.add_option("-e", "--strain", type="float",\
dest="strain",\
help="final strain to reach in simulation [%default]")
parser.add_option("--rate", type="float",\
dest="strainrate",\
help="(engineering) strain rate to simulate")
parser.add_option("-N", "--increments", type="int",\
dest="increments",\
help="number of increments to take")
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.add_option("-x", "--xmargin", type="float",\
dest="xmargin",\
help="margin in x in units of patch size [%default]")
parser.add_option("-y", "--ymargin", type="float",\
dest="ymargin",\
help="margin in y in units of patch size [%default]")
parser.add_option("-r", "--resolution", type="int",\
dest="resolution",\
help="number of Fourier points/Finite Elements across patch size + x_margin [%default]")
parser.add_option("-z", "--extrusion", type="int",\
dest="extrusion",\
help="number of repetitions in z-direction [%default]")
parser.add_option("-i", "--imagesize", type="int",\
dest="imgsize",\
help="size of PNG image")
parser.add_option("-M", "--coordtransformation", type="float", nargs=4, \
dest="M",\
help="2x2 transformation from rcb to Euler coords ( = M . [x_rcb,y_rcb])")
parser.add_option("--scatter", type="float",\
dest="scatter",\
help="orientation scatter [%default]")
parser.set_defaults(output = [])
parser.set_defaults(size = 1.0)
parser.set_defaults(xmargin = 0.0)
parser.set_defaults(ymargin = 0.0)
parser.set_defaults(resolution = 64)
parser.set_defaults(extrusion = 2)
parser.set_defaults(imgsize = 512)
parser.set_defaults(M = [0.0,1.0,1.0,0]) # M_11, M_12, M_21, M_22. x,y in RCB is y,x of Eulers!!
parser.set_defaults(tolerance = 1.0e-3)
parser.set_defaults(scatter = 0.0)
parser.set_defaults(strain = 0.2)
parser.set_defaults(strainrate = 1.0e-3)
parser.set_defaults(increments = 200)
parser.set_defaults(mesh = 'dt_planar_trimesh')
parser.set_defaults(twoD = False)
(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)
options.output = [s.lower() for s in options.output] # lower case
myName = os.path.splitext(args[0])[0]
print "\n%s\n"%myName
orientationData = rcbOrientationParser(boundarySegments)
rcData = rcbParser(boundarySegments,options.M,options.size,options.tolerance)
# ----- write image -----
if 'image' in options.output and options.imgsize > 0:
if ImageCapability:
image(myName,options.imgsize,options.xmargin,options.ymargin,rcData)
else:
print '...no image drawing possible (PIL missing)...'
# ----- write spectral geom -----
if 'spectral' in options.output:
fftdata = fftbuild(rcData, options.size, options.xmargin, options.ymargin, options.resolution, options.extrusion)
geomFile = open(myName+'_'+str(int(fftdata['resolution'][0]))+'.geom','w') # open geom file for writing
geomFile.write('resolution a %i b %i c %i\n'%(fftdata['resolution'])) # write resolution
geomFile.write('dimension x %f y %f z %f\n'%(fftdata['dimension'])) # write size
geomFile.write('homogenization 1\n') # write homogenization
geomFile.write('\n'.join(map(str,fftdata['fftpoints']))+'\n') # write grain indexes, one per line
geomFile.close() # close geom file
print('assigned %i out of %i Fourier points.'%(len(fftdata['fftpoints']), int(fftdata['resolution'][0])*int(fftdata['resolution'][1])))
# ----- write Mentat procedure -----
if 'mentat' in options.output:
if MentatCapability:
rcData['offsetPoints'] = 1+4 # gage definition generates 4 points
rcData['offsetSegments'] = 1+4 # gage definition generates 4 segments
cmds = [\
init(),
sample(options.size,rcData['dimension'][1]/rcData['dimension'][0],12,options.xmargin,options.ymargin),
patch(options.size,options.resolution,options.mesh,rcData),
gage(options.mesh,rcData),
]
if not options.twoD:
cmds += [expand3D(options.size*(1.0+2.0*options.xmargin)/options.resolution*options.extrusion,options.extrusion),]
cmds += [\
cleanUp(options.size),
materials(),
initial_conditions(len(rcData['grain']),rcData['grainMapping']),
boundary_conditions(options.strainrate,options.size*(1.0+2.0*options.xmargin)/options.resolution*options.extrusion,\
options.size,rcData['dimension'][1]/rcData['dimension'][0],options.xmargin,options.ymargin),
loadcase(options.strain/options.strainrate,options.increments,0.01),
job(len(rcData['grain']),rcData['grainMapping'],options.twoD),
postprocess(),
["*identify_sets","*regen","*fill_view","*save_as_model %s yes"%(myName)],
]
outputLocals = {'log':[]}
if (options.port != None):
py_connect('',options.port)
output(cmds,outputLocals,'Mentat')
py_disconnect()
if 'procedure' in options.output:
output(outputLocals['log'],outputLocals,'Stdout')
else:
print '...no interaction with Mentat possible...'
# ----- write config data to file -----
if 'mentat' in options.output or 'spectral' in options.output:
output = ''
output += '\n\n<homogenization>\n' + \
'\n[SX]\n' + \
'type\tisostrain\n' + \
'Ngrains\t1\n' + \
'\n\n<microstructure>\n'
for i,grain in enumerate(rcData['grainMapping']):
output += '\n[grain %i]\n'%grain + \
'crystallite\t1\n' + \
'(constituent)\tphase 1\ttexture %i\tfraction 1.0\n'%(i+1)
if (options.xmargin > 0.0):
output += '\n[x-margin]\n' + \
'crystallite\t1\n' + \
'(constituent)\tphase 2\ttexture %i\tfraction 1.0\n'%(len(rcData['grainMapping'])+1)
if (options.ymargin > 0.0):
output += '\n[y-margin]\n' + \
'crystallite\t1\n' + \
'(constituent)\tphase 2\ttexture %i\tfraction 1.0\n'%(len(rcData['grainMapping'])+1)
if (options.xmargin > 0.0 and options.ymargin > 0.0):
output += '\n[margin edge]\n' + \
'crystallite\t1\n' + \
'(constituent)\tphase 2\ttexture %i\tfraction 1.0\n'%(len(rcData['grainMapping'])+1)
output += '\n\n<crystallite>\n' + \
'\n[fillMeIn]\n' + \
'\n\n<phase>\n' + \
'\n[patch]\n'
if (options.xmargin > 0.0 or options.ymargin > 0.0):
output += '\n[margin]\n'
output += '\n\n<texture>\n\n'
for grain in rcData['grainMapping']:
output += '\n[grain %i]\n'%grain + \
'(gauss)\tphi1\t%f\tphi\t%f\tphi2\t%f\tscatter\t%f\tfraction\t1.0\n'\
%(math.degrees(orientationData[grain-1][0]),math.degrees(orientationData[grain-1][1]),math.degrees(orientationData[grain-1][2]),options.scatter)
if (options.xmargin > 0.0 or options.ymargin > 0.0):
output += '\n[margin]\n' + \
'(random)\t\tscatter\t0.0\tfraction\t1.0\n'
configFile = open(myName+'.config','w')
configFile.write(output)
configFile.close()