Operate on MFD file to apply servo links, not client-server communication

Speed up by many orders of magnitude... 64 x 64 x 64 mesh less than a minute! Thanks Brendan!

processing/pre/mentat_pbcOnBoxMesh.py

@@ -1,67 +1,103 @@
 #!/usr/bin/env python2.7
 # -*- coding: UTF-8 no BOM -*-
 
-import sys,os
+from __future__ import print_function
+import sys,os,re,time,tempfile
 import numpy as np
-from optparse import OptionParser
+import argparse
 import damask
 
+sys.path.append(damask.solver.Marc().libraryPath())
+
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 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([])
           
-active=[True,True,True] # directions on which to add PBC
+          # grab numbers
-def outMentat(cmd,locals):
+          nums = re.split(r'\s+', line.strip())
-  if cmd[0:3] == '(!)':
-    exec(cmd[3:])
-  elif cmd[0:3] == '(?)':
-    cmd = eval(cmd[3:])
-    py_mentat.py_send(cmd)
-  else:
-    py_mentat.py_send(cmd)
-  return
           
-#-------------------------------------------------------------------------------------------------
+          for num in nums:  
-def outFile(cmd,locals,dest):
+            # floating point has format ' -x.xxxxxxxxxxxxe+yy'
-  if cmd[0:3] == '(!)':
+            # scientific notation is used for float
-    exec(cmd[3:])
+            if (len(num) >= 4) and (num[-4] == 'e'):
-  elif cmd[0:3] == '(?)':
+              formatted[section]['els'][-1].append(float(num))
-    cmd = eval(cmd[3:])
+            else: # integer
-    dest.write(cmd+'\n')
+              formatted[section]['els'][-1].append(int(num))
-  else:
+        else: # not numerical data, so it is a label for an element or section end
-    dest.write(cmd+'\n')
+          if line[0] == '=' and re.search(r'=end=$', line) is not None: # End of section, avoiding regex if possible
-  return
+            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
-def output(cmds,locals,dest):
+
-  for cmd in cmds:
+def asMFD(mfd_data):
-    if isinstance(cmd,list):
+  result = ''
-      output(cmd,locals,dest)
+  for section in mfd_data:
-    else:
+    if section['uid'] > 0:
-      if dest == 'Mentat':
+      result += '=beg={0:5d} ({1})\n'.format(section['uid'], section['label'])
-        outMentat(str(cmd),locals)
+    for el in section['els']:
+      if type(el) == str:
+        result += el
+      elif type(el) == list:
+        for num in el:
+          if type(num) == int:
+            result += '{:20d}'.format(num)
+          elif type(num) == float:
+            result += '{:20.12e}'.format(num)
+          else:
+            damask.util.croak('WARNING: encountered unknown type: ' + str(type(el)))
+        result += '\n'
       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']
   box = {'min': np.zeros(3,dtype='d'),
          'max': 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['max'] = NodeCoords.max(axis=0)
   box['delta'] = box['max']-box['min']
@@ -108,9 +144,12 @@ def servoLink():
       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)]})
   
+  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
   
-  
+  links = {'uid': 1705, 'label': 'links', 'els': [[7,0],[9,0]]}
+  linkID = 0
   for node in linkNodes:                                                                            # loop over all linked nodes
     linkCoord = [node['coord']]                                                                     # start list of control node coords with my coords
     for dir in range(3):                                                                            # check for each direction
@@ -121,77 +160,100 @@ def servoLink():
 
     nLinks = len(linkCoord)
     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
 #-------------------------------------------------------------------------------------------------- 
-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)
 
-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)