Add MSC.Marc to VTK conversion support
marc_to_vtk.py takes Marc input file, subdivides it and creates a vtkUnstructuredGrid .vtu file with the resulting geometry. Currently supports hexahedron elements. vtk_addGridData.py adds nodal/cell data to some VTK grid. It is essentially vtk_addRectilinearGridData with support for unstructured grid (.vtu) format, which is marc_to_vtk.py's output.
This commit is contained in:
brendanvk
parent
beedd27a16
commit
a03bb1ec3e
|
|
@ -0,0 +1,167 @@
|
||||||
|
#!/usr/bin/env python2.7
|
||||||
|
# -*- coding: UTF-8 no BOM -*-
|
||||||
|
|
||||||
|
import os,re
|
||||||
|
import argparse
|
||||||
|
import damask
|
||||||
|
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.add_argument('filename', type=str, nargs='+', help='files to convert')
|
||||||
|
|
||||||
|
args = parser.parse_args()
|
||||||
|
files = args.filename
|
||||||
|
if type(files) is str:
|
||||||
|
files = [files]
|
||||||
|
|
||||||
|
|
||||||
|
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]
|
||||||
|
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
|
||||||
|
{
|
||||||
|
'type': int(line[10:20]),
|
||||||
|
'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]
|
||||||
|
coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL))
|
||||||
|
coordinates_header = coordinates_lines[0]
|
||||||
|
coordinates_lines = coordinates_lines[1:]
|
||||||
|
# marc input file does not use "e" in scientific notation, this adds it and converts
|
||||||
|
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([
|
||||||
|
fl_format(line[10:30]),
|
||||||
|
fl_format(line[30:50]),
|
||||||
|
fl_format(line[50:70])
|
||||||
|
])
|
||||||
|
), coordinates_lines))
|
||||||
|
|
||||||
|
# Subdivide volumes
|
||||||
|
grid = vtk.vtkUnstructuredGrid()
|
||||||
|
vertex_count = len(coordinates)
|
||||||
|
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
|
||||||
|
|
||||||
|
def subdivide_edge(vert1, vert2):
|
||||||
|
edge = ordered_pair(vert1, vert2)
|
||||||
|
|
||||||
|
if edge in edge_to_vert:
|
||||||
|
return edge_to_vert[edge]
|
||||||
|
|
||||||
|
newvert = len(coordinates) + 1
|
||||||
|
coordinates[newvert] = 0.5 * (coordinates[vert1] + coordinates[vert2]) # Average
|
||||||
|
edge_to_vert[edge] = newvert;
|
||||||
|
return newvert;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
for el_id in range(1, len(elements) + 1):
|
||||||
|
el = elements[el_id]
|
||||||
|
if el['type'] == 7:
|
||||||
|
# Hexahedron, subdivided
|
||||||
|
|
||||||
|
# There may be a better way to iterate over these, but this is consistent
|
||||||
|
# with the ordering scheme provided at https://damask.mpie.de/pub/Documentation/ElementType
|
||||||
|
|
||||||
|
subverts = np.zeros((3,3,3), dtype=int)
|
||||||
|
# Get corners
|
||||||
|
subverts[0, 0, 0] = el['verts'][0]
|
||||||
|
subverts[2, 0, 0] = el['verts'][1]
|
||||||
|
subverts[2, 2, 0] = el['verts'][2]
|
||||||
|
subverts[0, 2, 0] = el['verts'][3]
|
||||||
|
subverts[0, 0, 2] = el['verts'][4]
|
||||||
|
subverts[2, 0, 2] = el['verts'][5]
|
||||||
|
subverts[2, 2, 2] = el['verts'][6]
|
||||||
|
subverts[0, 2, 2] = el['verts'][7]
|
||||||
|
|
||||||
|
# lower edges
|
||||||
|
subverts[1, 0, 0] = subdivide_edge(subverts[0, 0, 0], subverts[2, 0, 0])
|
||||||
|
subverts[2, 1, 0] = subdivide_edge(subverts[2, 0, 0], subverts[2, 2, 0])
|
||||||
|
subverts[1, 2, 0] = subdivide_edge(subverts[2, 2, 0], subverts[0, 2, 0])
|
||||||
|
subverts[0, 1, 0] = subdivide_edge(subverts[0, 2, 0], subverts[0, 0, 0])
|
||||||
|
|
||||||
|
# middle edges
|
||||||
|
subverts[0, 0, 1] = subdivide_edge(subverts[0, 0, 0], subverts[0, 0, 2])
|
||||||
|
subverts[2, 0, 1] = subdivide_edge(subverts[2, 0, 0], subverts[2, 0, 2])
|
||||||
|
subverts[2, 2, 1] = subdivide_edge(subverts[2, 2, 0], subverts[2, 2, 2])
|
||||||
|
subverts[0, 2, 1] = subdivide_edge(subverts[0, 2, 0], subverts[0, 2, 2])
|
||||||
|
|
||||||
|
# top edges
|
||||||
|
subverts[1, 0, 2] = subdivide_edge(subverts[0, 0, 2], subverts[2, 0, 2])
|
||||||
|
subverts[2, 1, 2] = subdivide_edge(subverts[2, 0, 2], subverts[2, 2, 2])
|
||||||
|
subverts[1, 2, 2] = subdivide_edge(subverts[2, 2, 2], subverts[0, 2, 2])
|
||||||
|
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
|
||||||
|
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]
|
||||||
|
|
||||||
|
subverts[1, 0, 1] = subdivide_edge(subverts[1, 0, 0], subverts[1, 0, 2])
|
||||||
|
edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[2, 0, 1])] = subverts[1, 0, 1]
|
||||||
|
|
||||||
|
subverts[2, 1, 1] = subdivide_edge(subverts[2, 1, 0], subverts[2, 1, 2])
|
||||||
|
edge_to_vert[ordered_pair(subverts[2, 0, 1], subverts[2, 2, 1])] = subverts[2, 1, 1]
|
||||||
|
|
||||||
|
subverts[1, 2, 1] = subdivide_edge(subverts[1, 2, 0], subverts[1, 2, 2])
|
||||||
|
edge_to_vert[ordered_pair(subverts[0, 2, 1], subverts[2, 2, 1])] = subverts[1, 2, 1]
|
||||||
|
|
||||||
|
subverts[0, 1, 1] = subdivide_edge(subverts[0, 1, 0], subverts[0, 1, 2])
|
||||||
|
edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[0, 2, 1])] = subverts[0, 1, 1]
|
||||||
|
|
||||||
|
subverts[1, 1, 2] = subdivide_edge(subverts[1, 0, 2], subverts[1, 2, 2])
|
||||||
|
edge_to_vert[ordered_pair(subverts[0, 1, 2], subverts[2, 1, 2])] = subverts[1, 1, 2]
|
||||||
|
|
||||||
|
# and finally the center. There are three ways to calculate, but elements should
|
||||||
|
# not intersect, so the edge_to_vert part isn't needed here.
|
||||||
|
subverts[1, 1, 1] = subdivide_edge(subverts[1, 1, 0], subverts[1, 1, 2])
|
||||||
|
|
||||||
|
|
||||||
|
# Now make the hexahedron subelements
|
||||||
|
# order in which vtk expects vertices for a hexahedron
|
||||||
|
order = np.array([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)])
|
||||||
|
for z in range(2):
|
||||||
|
for y in range(2):
|
||||||
|
for x in range(2):
|
||||||
|
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
|
||||||
|
grid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds())
|
||||||
|
|
||||||
|
|
||||||
|
else:
|
||||||
|
damask.util.croak('Unsupported Marc element type: {} (skipping)'.format(el['type']))
|
||||||
|
|
||||||
|
# Load all points
|
||||||
|
points = vtk.vtkPoints()
|
||||||
|
for point in range(1, len(coordinates) + 1): # marc indices start at 1
|
||||||
|
points.InsertNextPoint(coordinates[point].tolist())
|
||||||
|
|
||||||
|
grid.SetPoints(points)
|
||||||
|
|
||||||
|
# grid now contains the elements from the given marc file
|
||||||
|
writer = vtk.vtkXMLUnstructuredGridWriter()
|
||||||
|
writer.SetFileName(re.sub(r'\..+', ".vtu", f)) # *.vtk extension does not work in paraview
|
||||||
|
#writer.SetCompressorTypeToZLib()
|
||||||
|
|
||||||
|
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
|
||||||
|
else: writer.SetInputData(grid)
|
||||||
|
writer.Write()
|
||||||
|
|
@ -0,0 +1,206 @@
|
||||||
|
#!/usr/bin/env python2.7
|
||||||
|
# -*- coding: UTF-8 no BOM -*-
|
||||||
|
|
||||||
|
import os,vtk
|
||||||
|
import damask
|
||||||
|
from vtk.util import numpy_support
|
||||||
|
from collections import defaultdict
|
||||||
|
from optparse import OptionParser
|
||||||
|
|
||||||
|
scriptName = os.path.splitext(os.path.basename(__file__))[0]
|
||||||
|
scriptID = ' '.join([scriptName,damask.version])
|
||||||
|
|
||||||
|
# --------------------------------------------------------------------
|
||||||
|
# MAIN
|
||||||
|
# --------------------------------------------------------------------
|
||||||
|
|
||||||
|
msg = "Add scalars, vectors, and/or an RGB tuple from"
|
||||||
|
msg += "an ASCIItable to existing VTK grid (.vtr/.vtk/.vtu)."
|
||||||
|
parser = OptionParser(option_class=damask.extendableOption,
|
||||||
|
usage='%prog options [file[s]]',
|
||||||
|
description = msg,
|
||||||
|
version = scriptID)
|
||||||
|
|
||||||
|
parser.add_option( '--vtk',
|
||||||
|
dest = 'vtk',
|
||||||
|
type = 'string', metavar = 'string',
|
||||||
|
help = 'VTK file name')
|
||||||
|
parser.add_option( '--inplace',
|
||||||
|
dest = 'inplace',
|
||||||
|
action = 'store_true',
|
||||||
|
help = 'modify VTK file in-place')
|
||||||
|
parser.add_option('-r', '--render',
|
||||||
|
dest = 'render',
|
||||||
|
action = 'store_true',
|
||||||
|
help = 'open output in VTK render window')
|
||||||
|
parser.add_option('-d', '--data',
|
||||||
|
dest = 'data',
|
||||||
|
action = 'extend', metavar = '<string LIST>',
|
||||||
|
help = 'scalar/vector value(s) label(s)')
|
||||||
|
parser.add_option('-t', '--tensor',
|
||||||
|
dest = 'tensor',
|
||||||
|
action = 'extend', metavar = '<string LIST>',
|
||||||
|
help = 'tensor (3x3) value label(s)')
|
||||||
|
parser.add_option('-c', '--color',
|
||||||
|
dest = 'color',
|
||||||
|
action = 'extend', metavar = '<string LIST>',
|
||||||
|
help = 'RGB color tuple label')
|
||||||
|
|
||||||
|
parser.set_defaults(data = [],
|
||||||
|
tensor = [],
|
||||||
|
color = [],
|
||||||
|
inplace = False,
|
||||||
|
render = False,
|
||||||
|
)
|
||||||
|
|
||||||
|
(options, filenames) = parser.parse_args()
|
||||||
|
|
||||||
|
if not options.vtk: parser.error('No VTK file specified.')
|
||||||
|
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
|
||||||
|
|
||||||
|
if os.path.splitext(options.vtk)[1] == '.vtr':
|
||||||
|
reader = vtk.vtkXMLRectilinearGridReader()
|
||||||
|
reader.SetFileName(options.vtk)
|
||||||
|
reader.Update()
|
||||||
|
rGrid = reader.GetOutput()
|
||||||
|
writer = vtk.vtkXMLRectilinearGridWriter()
|
||||||
|
writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr'))
|
||||||
|
elif os.path.splitext(options.vtk)[1] == '.vtk':
|
||||||
|
reader = vtk.vtkGenericDataObjectReader()
|
||||||
|
reader.SetFileName(options.vtk)
|
||||||
|
reader.Update()
|
||||||
|
rGrid = reader.GetRectilinearGridOutput()
|
||||||
|
writer = vtk.vtkXMLRectilinearGridWriter()
|
||||||
|
writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr'))
|
||||||
|
elif os.path.splitext(options.vtk)[1] == '.vtu':
|
||||||
|
reader = vtk.vtkXMLUnstructuredGridReader()
|
||||||
|
reader.SetFileName(options.vtk)
|
||||||
|
reader.Update()
|
||||||
|
rGrid = reader.GetOutput()
|
||||||
|
writer = vtk.vtkXMLUnstructuredGridWriter()
|
||||||
|
writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtu' if options.inplace else '_added.vtu'))
|
||||||
|
else:
|
||||||
|
parser.error('Unsupported VTK file type extension.')
|
||||||
|
|
||||||
|
Npoints = rGrid.GetNumberOfPoints()
|
||||||
|
Ncells = rGrid.GetNumberOfCells()
|
||||||
|
|
||||||
|
damask.util.croak('{}: {} points and {} cells...'.format(options.vtk,Npoints,Ncells))
|
||||||
|
|
||||||
|
# --- loop over input files -------------------------------------------------------------------------
|
||||||
|
|
||||||
|
if filenames == []: filenames = [None]
|
||||||
|
|
||||||
|
for name in filenames:
|
||||||
|
try: table = damask.ASCIItable(name = name,
|
||||||
|
buffered = False,
|
||||||
|
readonly = True)
|
||||||
|
except: continue
|
||||||
|
damask.util.report(scriptName, name)
|
||||||
|
|
||||||
|
# --- interpret header ----------------------------------------------------------------------------
|
||||||
|
|
||||||
|
table.head_read()
|
||||||
|
|
||||||
|
remarks = []
|
||||||
|
errors = []
|
||||||
|
VTKarray = {}
|
||||||
|
active = defaultdict(list)
|
||||||
|
|
||||||
|
for datatype,dimension,label in [['data',99,options.data],
|
||||||
|
['tensor',9,options.tensor],
|
||||||
|
['color' ,3,options.color],
|
||||||
|
]:
|
||||||
|
for i,dim in enumerate(table.label_dimension(label)):
|
||||||
|
me = label[i]
|
||||||
|
if dim == -1: remarks.append('{} "{}" not found...'.format(datatype,me))
|
||||||
|
elif dim > dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension))
|
||||||
|
else:
|
||||||
|
remarks.append('adding {} "{}"...'.format(datatype,me))
|
||||||
|
active[datatype].append(me)
|
||||||
|
|
||||||
|
if remarks != []: damask.util.croak(remarks)
|
||||||
|
if errors != []:
|
||||||
|
damask.util.croak(errors)
|
||||||
|
table.close(dismiss = True)
|
||||||
|
continue
|
||||||
|
|
||||||
|
# ------------------------------------------ process data ---------------------------------------
|
||||||
|
|
||||||
|
table.data_readArray([item for sublist in active.values() for item in sublist]) # read all requested data
|
||||||
|
|
||||||
|
for datatype,labels in active.items(): # loop over scalar,color
|
||||||
|
for me in labels: # loop over all requested items
|
||||||
|
VTKtype = vtk.VTK_DOUBLE
|
||||||
|
VTKdata = table.data[:, table.label_indexrange(me)].copy() # copy to force contiguous layout
|
||||||
|
|
||||||
|
if datatype == 'color':
|
||||||
|
VTKtype = vtk.VTK_UNSIGNED_CHAR
|
||||||
|
VTKdata = (VTKdata*255).astype(int) # translate to 0..255 UCHAR
|
||||||
|
elif datatype == 'tensor':
|
||||||
|
VTKdata[:,1] = VTKdata[:,3] = 0.5*(VTKdata[:,1]+VTKdata[:,3])
|
||||||
|
VTKdata[:,2] = VTKdata[:,6] = 0.5*(VTKdata[:,2]+VTKdata[:,6])
|
||||||
|
VTKdata[:,5] = VTKdata[:,7] = 0.5*(VTKdata[:,5]+VTKdata[:,7])
|
||||||
|
|
||||||
|
VTKarray[me] = numpy_support.numpy_to_vtk(num_array=VTKdata,deep=True,array_type=VTKtype)
|
||||||
|
VTKarray[me].SetName(me)
|
||||||
|
|
||||||
|
table.close() # close input ASCII table
|
||||||
|
|
||||||
|
# ------------------------------------------ add data ---------------------------------------
|
||||||
|
|
||||||
|
if len(table.data) == Npoints: mode = 'point'
|
||||||
|
elif len(table.data) == Ncells: mode = 'cell'
|
||||||
|
else:
|
||||||
|
damask.util.croak('Data count is incompatible with grid...')
|
||||||
|
continue
|
||||||
|
|
||||||
|
damask.util.croak('{} mode...'.format(mode))
|
||||||
|
|
||||||
|
for datatype,labels in active.items(): # loop over scalar,color
|
||||||
|
if datatype == 'color':
|
||||||
|
if mode == 'cell': rGrid.GetCellData().SetScalars(VTKarray[active['color'][0]])
|
||||||
|
elif mode == 'point': rGrid.GetPointData().SetScalars(VTKarray[active['color'][0]])
|
||||||
|
for me in labels: # loop over all requested items
|
||||||
|
if mode == 'cell': rGrid.GetCellData().AddArray(VTKarray[me])
|
||||||
|
elif mode == 'point': rGrid.GetPointData().AddArray(VTKarray[me])
|
||||||
|
|
||||||
|
rGrid.Modified()
|
||||||
|
if vtk.VTK_MAJOR_VERSION <= 5: rGrid.Update()
|
||||||
|
|
||||||
|
# ------------------------------------------ output result ---------------------------------------
|
||||||
|
|
||||||
|
writer.SetDataModeToBinary()
|
||||||
|
writer.SetCompressorTypeToZLib()
|
||||||
|
if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
|
||||||
|
else: writer.SetInputData(rGrid)
|
||||||
|
writer.Write()
|
||||||
|
|
||||||
|
# ------------------------------------------ render result ---------------------------------------
|
||||||
|
|
||||||
|
if options.render:
|
||||||
|
mapper = vtk.vtkDataSetMapper()
|
||||||
|
mapper.SetInputData(rGrid)
|
||||||
|
actor = vtk.vtkActor()
|
||||||
|
actor.SetMapper(mapper)
|
||||||
|
|
||||||
|
# Create the graphics structure. The renderer renders into the
|
||||||
|
# render window. The render window interactor captures mouse events
|
||||||
|
# and will perform appropriate camera or actor manipulation
|
||||||
|
# depending on the nature of the events.
|
||||||
|
|
||||||
|
ren = vtk.vtkRenderer()
|
||||||
|
|
||||||
|
renWin = vtk.vtkRenderWindow()
|
||||||
|
renWin.AddRenderer(ren)
|
||||||
|
|
||||||
|
ren.AddActor(actor)
|
||||||
|
ren.SetBackground(1, 1, 1)
|
||||||
|
renWin.SetSize(200, 200)
|
||||||
|
|
||||||
|
iren = vtk.vtkRenderWindowInteractor()
|
||||||
|
iren.SetRenderWindow(renWin)
|
||||||
|
|
||||||
|
iren.Initialize()
|
||||||
|
renWin.Render()
|
||||||
|
iren.Start()
|
||||||
Loading…
Reference in New Issue