Subdivide vertex data with geometry when converting Marc input to VTK

processing/post/marc_to_vtk.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python2.7
 # -*- coding: UTF-8 no BOM -*-
 
-import os,re
+import os,sys,re
 import argparse
 import damask
 import vtk, numpy as np
@@ -9,25 +9,39 @@ import vtk, numpy as np
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName, damask.version])
 
-parser = argparse.ArgumentParser(description='Convert from Marc input file format to VTK', version = scriptID)
+parser = argparse.ArgumentParser(description='Convert from Marc input file format (.dat) to VTK format (.vtu)', version = scriptID)
-parser.add_argument('filename', type=str, nargs='+', help='files to convert')
+parser.add_argument('filename', type=str, help='file to convert')
+parser.add_argument('-t', '--table',  type=str, help='ASCIItable file containing nodal data to subdivide and interpolate')
 
 args = parser.parse_args()
-files = args.filename
-if type(files) is str:
-  files = [files]
 
-    
+with open(args.filename, 'r') as marcfile:
-for f in files:
-  with open(f, 'r') as marcfile:
   marctext = marcfile.read();
-  # Extract connectivity chunk from file...
+  
-  connectivity_text = re.findall(r'connectivity[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+# Load table (if any)
-  connectivity_lines = re.split(r'[\n\r]+', connectivity_text, flags=(re.MULTILINE | re.DOTALL))
+if args.table is not None:
-  connectivity_header = connectivity_lines[0]
+  try:
-  connectivity_lines = connectivity_lines[1:]
+    table = damask.ASCIItable(
-  # Construct element map
+      name=args.table,
-  elements = dict(map(lambda line: 
+      outname='subdivided_{}'.format(args.table),
+      buffered=True
+    )
+    
+    table.head_read()
+    table.data_readArray()
+    
+    # Python list is faster for appending
+    nodal_data = list(table.data)
+  except: args.table = None
+
+# Extract connectivity chunk from file...
+connectivity_text = re.findall(r'connectivity[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+connectivity_lines = re.split(r'[\n\r]+', connectivity_text, flags=(re.MULTILINE | re.DOTALL))
+connectivity_header = connectivity_lines[0]
+connectivity_lines = connectivity_lines[1:]
+
+# Construct element map
+elements = dict(map(lambda line: 
   (
     int(line[0:10]), # index
     { 
@@ -35,15 +49,17 @@ for f in files:
       'verts': list(map(int, re.split(r' +', line[20:].strip())))
     }
   ), connectivity_lines))
-  # Extract coordinate chunk from file
+
-  coordinates_text = re.findall(r'coordinates[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+# Extract coordinate chunk from file
-  coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL))
+coordinates_text = re.findall(r'coordinates[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
-  coordinates_header = coordinates_lines[0]
+coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL))
-  coordinates_lines = coordinates_lines[1:]
+coordinates_header = coordinates_lines[0]
-  # marc input file does not use "e" in scientific notation, this adds it and converts
+coordinates_lines = coordinates_lines[1:]
-  fl_format = lambda string: float(re.sub(r'(\d)([\+\-])', r'\1e\2', string))
+
-  # Construct coordinate map
+# marc input file does not use "e" in scientific notation, this adds it and converts
-  coordinates = dict(map(lambda line:
+fl_format = lambda string: float(re.sub(r'(\d)([\+\-])', r'\1e\2', string))
+# Construct coordinate map
+coordinates = dict(map(lambda line:
   (
     int(line[0:10]),
     np.array([
@@ -53,26 +69,29 @@ for f in files:
     ])
   ), coordinates_lines))
 
-  # Subdivide volumes
+# Subdivide volumes
-  grid = vtk.vtkUnstructuredGrid()
+grid = vtk.vtkUnstructuredGrid()
-  vertex_count = len(coordinates)
+vertex_count = len(coordinates)
-  edge_to_vert = dict() # when edges are subdivided, a new vertex in the middle is produced and placed in here
+edge_to_vert = dict() # when edges are subdivided, a new vertex in the middle is produced and placed in here
-  ordered_pair = lambda a, b: (a, b) if a < b else (b, a) # edges are bidirectional
+ordered_pair = lambda a, b: (a, b) if a < b else (b, a) # edges are bidirectional
 
-  def subdivide_edge(vert1, vert2):
+def subdivide_edge(vert1, vert2):
   edge = ordered_pair(vert1, vert2)
   
   if edge in edge_to_vert:
     return edge_to_vert[edge]
   
+  # Vertex does not exist, create it
   newvert = len(coordinates) + 1
   coordinates[newvert] = 0.5 * (coordinates[vert1] + coordinates[vert2]) # Average
   edge_to_vert[edge] = newvert;
+  
+  # Interpolate nodal data
+  if args.table is not None:
+    nodal_data.append(0.5 * (nodal_data[vert1 - 1] + nodal_data[vert2 - 1]))
   return newvert;
 
-  
+for el_id in range(1, len(elements) + 1): # Marc starts counting at 1
-  
-  for el_id in range(1, len(elements) + 1):
   el = elements[el_id]
   if el['type'] == 7:
     # Hexahedron, subdivided
@@ -110,7 +129,8 @@ for f in files:
     subverts[0, 1, 2] = subdivide_edge(subverts[0, 2, 2], subverts[0, 0, 2])
     
     # then faces...  The edge_to_vert addition is due to there being two ways
-      # to calculate a face, depending which opposite vertices are used to subdivide
+    # to calculate a face vertex, depending on which opposite vertices are used to subdivide.
+    # This way, we avoid creating duplicate vertices.
     subverts[1, 1, 0] = subdivide_edge(subverts[1, 0, 0], subverts[1, 2, 0])
     edge_to_vert[ordered_pair(subverts[0, 1, 0], subverts[2, 1, 0])] = subverts[1, 1, 0]
     
@@ -143,25 +163,37 @@ for f in files:
           hex_ = vtk.vtkHexahedron()
           for vert_id in range(8):
             coord = order[vert_id] + (x, y, z)
-              hex_.GetPointIds().SetId(vert_id, subverts[coord[0], coord[1], coord[2]] - 1) # minus one, since vtk starts at zero but marc starts at one
+            # minus one, since vtk starts at zero but marc starts at one
+            hex_.GetPointIds().SetId(vert_id, subverts[coord[0], coord[1], coord[2]] - 1) 
           grid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds())
     
-      
   else:
     damask.util.croak('Unsupported Marc element type: {} (skipping)'.format(el['type']))
 
-  # Load all points
+# Load all points
-  points = vtk.vtkPoints()
+points = vtk.vtkPoints()
-  for point in range(1, len(coordinates) + 1): # marc indices start at 1
+for point in range(1, len(coordinates) + 1): # marc indices start at 1
   points.InsertNextPoint(coordinates[point].tolist())
 
-  grid.SetPoints(points)
+grid.SetPoints(points)
 
-  # grid now contains the elements from the given marc file
+# grid now contains the elements from the given marc file
-  writer = vtk.vtkXMLUnstructuredGridWriter()
+writer = vtk.vtkXMLUnstructuredGridWriter()
-  writer.SetFileName(re.sub(r'\..+', ".vtu", f))  # *.vtk extension does not work in paraview
+writer.SetFileName(re.sub(r'\..+', ".vtu", args.filename))  # *.vtk extension does not work in paraview
-  #writer.SetCompressorTypeToZLib()
 
-  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
+if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
-  else:                          writer.SetInputData(grid)
+else:                          writer.SetInputData(grid)
-  writer.Write()
+writer.Write()
+
+if args.table is not None:  
+  table.info_append([
+    scriptID + ' ' + ' '.join(sys.argv[1:]),
+  ])
+  table.head_write()
+  table.output_flush()
+  
+  table.data = np.array(nodal_data)
+  
+  table.data_writeArray()
+
+table.close()