Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Franz Roters 2017-04-27 08:22:08 +02:00
commit dd34f28cdb
2 changed files with 171 additions and 109 deletions

View File

@ -1 +1 @@
v2.0.1-626-gc10c58b v2.0.1-631-g175e119

278
processing/pre/mentat_pbcOnBoxMesh.py Executable file → Normal file
View File

@ -1,67 +1,103 @@
#!/usr/bin/env python2.7 #!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import sys,os from __future__ import print_function
import sys,os,re,time,tempfile
import numpy as np import numpy as np
from optparse import OptionParser import argparse
import damask import damask
sys.path.append(damask.solver.Marc().libraryPath())
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
sys.path.append(damask.solver.Marc().libraryPath()) # Convert .mfd file into a usable format
# Broken into labeled sections (eg. nodes, links, etc)
# Each section has a list of labeled elements with formatted numerical data
def parseMFD(fname):
formatted = []
with open(fname,'r') as dat:
section = 0
formatted.append({'label': 'header', 'uid': -1, 'els': []})
# in between =beg= and =end= part of file
in_block = False
for line in dat:
if in_block: # currently in a section
# lines that start with a space are numerical data
if line[0] == ' ':
formatted[section]['els'].append([])
# grab numbers
nums = re.split(r'\s+', line.strip())
for num in nums:
# floating point has format ' -x.xxxxxxxxxxxxe+yy'
# scientific notation is used for float
if (len(num) >= 4) and (num[-4] == 'e'):
formatted[section]['els'][-1].append(float(num))
else: # integer
formatted[section]['els'][-1].append(int(num))
else: # not numerical data, so it is a label for an element or section end
if line[0] == '=' and re.search(r'=end=$', line) is not None: # End of section, avoiding regex if possible
in_block = False
else:
formatted[section]['els'].append([])
formatted[section]['els'][-1] = line
else: # Not in a section, we are looking for a =beg= now
search = re.search(r'=beg=\s+(\d+)\s\((.*?)\)', line)
if search is not None: # found start of a new section
section += 1
in_block = True
formatted.append({'label': search.group(2), 'uid': int(search.group(1)), 'els': []})
else: # No =beg= found, probably in the header
# Either header or somthing we didn't plan for - just save the line so it isn't lost
if formatted[section]['uid'] > 0:
section += 1
formatted.append({'label': '', 'uid': -2, 'els': []}) # make dummy section to store unrecognized data
formatted[section]['els'].append(line)
return formatted
active=[True,True,True] # directions on which to add PBC def asMFD(mfd_data):
def outMentat(cmd,locals): result = ''
if cmd[0:3] == '(!)': for section in mfd_data:
exec(cmd[3:]) if section['uid'] > 0:
elif cmd[0:3] == '(?)': result += '=beg={0:5d} ({1})\n'.format(section['uid'], section['label'])
cmd = eval(cmd[3:]) for el in section['els']:
py_mentat.py_send(cmd) if type(el) == str:
else: result += el
py_mentat.py_send(cmd) elif type(el) == list:
return for num in el:
if type(num) == int:
#------------------------------------------------------------------------------------------------- result += '{:20d}'.format(num)
def outFile(cmd,locals,dest): elif type(num) == float:
if cmd[0:3] == '(!)': result += '{:20.12e}'.format(num)
exec(cmd[3:]) else:
elif cmd[0:3] == '(?)': damask.util.croak('WARNING: encountered unknown type: ' + str(type(el)))
cmd = eval(cmd[3:]) result += '\n'
dest.write(cmd+'\n')
else:
dest.write(cmd+'\n')
return
#-------------------------------------------------------------------------------------------------
def output(cmds,locals,dest):
for cmd in cmds:
if isinstance(cmd,list):
output(cmd,locals,dest)
else:
if dest == 'Mentat':
outMentat(str(cmd),locals)
else: else:
outFile(str(cmd),locals,dest) damask.util.croak('WARNING: encountered unknown type: ' + str(type(el)))
return if section['uid'] > 0:
result += '=end=\n'
return result.strip()
def add_servoLinks(mfd_data,active=[True,True,True]): # directions on which to add PBC
def servoLink():
cmds = []
base = ['x','y','z'] base = ['x','y','z']
box = {'min': np.zeros(3,dtype='d'), box = {'min': np.zeros(3,dtype='d'),
'max': np.zeros(3,dtype='d'), 'max': np.zeros(3,dtype='d'),
'delta': np.zeros(3,dtype='d'), 'delta': np.zeros(3,dtype='d'),
} }
Nnodes = py_mentat.py_get_int("nnodes()")
NodeCoords = np.zeros((Nnodes,3),dtype='d') mfd_dict = {}
for node in range(Nnodes): for i in range(len(mfd_data)):
NodeCoords[node,0] = py_mentat.py_get_float("node_x(%i)"%(node+1)) mfd_dict[mfd_data[i]['label']] = i
NodeCoords[node,1] = py_mentat.py_get_float("node_y(%i)"%(node+1))
NodeCoords[node,2] = py_mentat.py_get_float("node_z(%i)"%(node+1)) NodeCoords = np.array(mfd_data[mfd_dict['nodes']]['els'][0::4])[:,1:4]
Nnodes = NodeCoords.shape[0]
box['min'] = NodeCoords.min(axis=0) # find the bounding box box['min'] = NodeCoords.min(axis=0) # find the bounding box
box['max'] = NodeCoords.max(axis=0) box['max'] = NodeCoords.max(axis=0)
box['delta'] = box['max']-box['min'] box['delta'] = box['max']-box['min']
@ -108,9 +144,12 @@ def servoLink():
if any([maxFlag[i] and active[i] for i in range(3)]): if any([maxFlag[i] and active[i] for i in range(3)]):
linkNodes.append({'id': node+1,'coord': NodeCoords[node], 'faceMember': [maxFlag[i] and active[i] for i in range(3)]}) linkNodes.append({'id': node+1,'coord': NodeCoords[node], 'faceMember': [maxFlag[i] and active[i] for i in range(3)]})
mfd_data[mfd_dict['entities']]['els'][0][0] += len(linkNodes) * 3
baseCorner = baseNode["%.8e"%box['min'][0]]["%.8e"%box['min'][1]]["%.8e"%box['min'][2]] # detect ultimate base node baseCorner = baseNode["%.8e"%box['min'][0]]["%.8e"%box['min'][1]]["%.8e"%box['min'][2]] # detect ultimate base node
links = {'uid': 1705, 'label': 'links', 'els': [[7,0],[9,0]]}
linkID = 0
for node in linkNodes: # loop over all linked nodes for node in linkNodes: # loop over all linked nodes
linkCoord = [node['coord']] # start list of control node coords with my coords linkCoord = [node['coord']] # start list of control node coords with my coords
for dir in range(3): # check for each direction for dir in range(3): # check for each direction
@ -121,77 +160,100 @@ def servoLink():
nLinks = len(linkCoord) nLinks = len(linkCoord)
for dof in [1,2,3]: for dof in [1,2,3]:
cmds.append([ tied_node = node['id']
"*new_link *link_class servo", nterms = 1 + nLinks
"*link_class servo *tied_node %i"%node['id'],
"*link_class servo *tied_dof %i"%dof, linkID += 1
"*servo_nterms %i"%(1+nLinks), # Link header
]) links['els'].append('link{0}\n'.format(linkID))
for i in range(nLinks): links['els'].append([linkID, 1])
cmds.append([ links['els'].append([0])
"*link_class servo *servo_ret_node %i %i"\ links['els'].append([0])
%(i+1,baseNode["%.8e"%linkCoord[i][0]]["%.8e"%linkCoord[i][1]]["%.8e"%linkCoord[i][2]]), links['els'].append([0, 0, 0, tied_node])
"*link_class servo *servo_ret_dof %i %i"%(i+1,dof),
"*link_class servo *servo_ret_coef %i 1"%(i+1), # these need to be put in groups of four
]) link_payload = [dof, 0, nterms]
cmds.append([
"*link_class servo *servo_ret_node %i %i"%(1+nLinks,baseCorner), # Individual node contributions (node, dof, coef.)
"*link_class servo *servo_ret_dof %i %i"%(1+nLinks,dof), for i in range(nterms):
"*link_class servo *servo_ret_coef %i -%i"%(1+nLinks,nLinks-1), if i == nLinks:
]) link_payload.append(baseCorner)
return cmds else:
link_payload.append(baseNode["%.8e"%linkCoord[i][0]]["%.8e"%linkCoord[i][1]]["%.8e"%linkCoord[i][2]])
for i in range(nterms):
link_payload.append(dof)
for i in range(nterms):
if i == nLinks:
link_payload.append(1.0 - nLinks)
else:
link_payload.append(1.0)
# Needs to be formatted 4 data points per row, character width of 20, so 80 total
for j in range(0, len(link_payload), 4):
links['els'].append(link_payload[j:j+4])
if j+4 < len(link_payload):
links['els'].append(link_payload[j+4:])
i = 0
while i < len(mfd_data) and mfd_data[i]['uid'] < 1705: i += 1
if mfd_data[i]['uid'] == 1705: del mfd_data[i]
mfd_data.insert(i, links)
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# MAIN # MAIN
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage = '%prog [options]', description = """ parser = argparse.ArgumentParser(description = """
Set up servo linking to achieve periodic boundary conditions for a regular hexahedral mesh presently opened in MSC.Mentat Set up servo linking to achieve periodic boundary conditions for a regular hexahedral mesh.
Use *py_connection to operate on model presently opened in MSC.Mentat.
""", version = scriptID) """, version = scriptID)
parser.add_option("-p", "--port", type="int", dest="port", metavar='int', parser.add_argument('-p', '--port',
help="Mentat connection port [%default]") type = int, metavar = 'int', default = None,
parser.add_option("-v", "--verbose", action="store_true", dest="verbose", help = 'Mentat connection port')
help="write Mentat command stream also to stdout [%default]") parser.add_argument('-x',
parser.set_defaults(port = 40007, action = 'store_false',
verbose = False) help = 'no PBC along x direction')
parser.add_argument('-y',
action = 'store_false',
help = 'no PBC along y direction')
parser.add_argument('-z',
action = 'store_false',
help = 'no PBC along z direction')
parser.add_argument('file', nargs='*',
help = 'Mentat formatted data (.mfd) file[s] to add periodic boundary conditions to')
(options, args) = parser.parse_args() args = parser.parse_args()
try: remote = args.port is not None
import py_mentat
except:
damask.util.croak('error: no valid Mentat release found')
sys.exit(-1)
outputLocals = {} if remote:
try: import py_mentat
except:
damask.util.croak('no valid Mentat release found.')
sys.exit(-1)
damask.util.report(scriptName,'waiting to connect...') damask.util.report(scriptName, 'waiting to connect...')
args.file = [os.path.join(tempfile._get_default_tempdir(), next(tempfile._get_candidate_names()) + '.mfd')]
try:
py_mentat.py_connect('',args.port)
py_mentat.py_send('*set_save_formatted on')
py_mentat.py_send('*save_as_model "{}" yes'.format(args.file[0]))
py_mentat.py_get_int("nnodes()") # hopefully blocks until file is written
except:
damask.util.croak('failed. try setting Tools/Python/"Run as Separate Process" & "Initiate".')
sys.exit()
damask.util.croak( 'connected...')
try: for mfdfile in args.file:
py_mentat.py_connect('',options.port) while remote and not os.path.exists(mfdfile): time.sleep(0.5) # wait for Mentat to write MFD file
# prevent redrawing in Mentat, should be much faster. damask.util.report(scriptName, mfdfile)
# Since py_connect has no return value, try this to determine if failed or not mfd = parseMFD(mfdfile)
output(['*draw_manual'],outputLocals,'Mentat') add_servoLinks(mfd,[args.x,args.y,args.z])
except: with open(mfdfile, 'w') as file:
damask.util.croak('Could not connect. Set Tools/Python/"Run as Separate Process" & "Initiate"...') file.write(asMFD(mfd))
sys.exit()
damask.util.croak('connected...') if remote:
try: py_mentat.py_send('*open_model "{}"'.format(args.file[0]))
output(['*remove_all_servos', except: damask.util.croak('lost connection on sending open command for "{}".'.format(args.file[0]))
'*sweep_all',
'*renumber_nodes',
'*set_links off',
],outputLocals,'Mentat') # script depends on consecutive numbering of nodes
cmds = servoLink()
output(cmds,outputLocals,'Mentat')
output(['*draw_automatic',
],outputLocals,'Mentat') # script depends on consecutive numbering of nodes
py_mentat.py_disconnect()
if options.verbose:
output(cmds,outputLocals,sys.stdout)