223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

vtk files, mainly removed from vtk import *

This commit is contained in:
Martin Diehl 2016-03-02 13:52:33 +01:00
parent f7fedc4744
commit 9fa49b8584
3 changed files with 65 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

118
processing/post/vtk2ang.py
View File

@ -4,17 +4,15 @@
import os,string,math,sys import os,string,math,sys
import numpy as np import numpy as np
from optparse import OptionParser from optparse import OptionParser
from vtk import * import vtk
import damask
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])
# ----------------------------- # -----------------------------
def getHeader(filename,sizeFastIndex,sizeSlowIndex,stepsize): def getHeader(filename,sizeFastIndex,sizeSlowIndex,stepsize):
# ----------------------------- """returns header for ang file step size in micrometer"""
# returns header for ang file
# step size in micrometer
return '\n'.join([ \ return '\n'.join([ \
'# TEM_PIXperUM 1.000000', \ '# TEM_PIXperUM 1.000000', \
'# x-star 1.000000', \ '# x-star 1.000000', \
@ -50,10 +48,7 @@ def getHeader(filename,sizeFastIndex,sizeSlowIndex,stepsize):
# ----------------------------- # -----------------------------
def positiveRadians(angle): def positiveRadians(angle):
# ----------------------------- """returns positive angle in radians from angle in degrees"""
# returns positive angle in radians
# gets angle in degrees
angle = math.radians(float(angle)) angle = math.radians(float(angle))
while angle < 0.0: while angle < 0.0:
angle += 2.0 * math.pi angle += 2.0 * math.pi
@ -63,11 +58,12 @@ def positiveRadians(angle):
# ----------------------------- # -----------------------------
def getDataLine(angles,x,y,validData=True): def getDataLine(angles,x,y,validData=True):
# ----------------------------- """
# returns string of one line in ang file returns string of one line in ang file
# convention in ang file: y coordinate comes first and is fastest index
# positions in micrometer convention in ang file: y coordinate comes first and is fastest index
positions in micrometer
"""
info = {True: (9999.9, 1.0, 0,99999,0.0), info = {True: (9999.9, 1.0, 0,99999,0.0),
False: ( -1.0,-1.0,-1, -1,1.0)} False: ( -1.0,-1.0,-1, -1,1.0)}
return '%9.5f %9.5f %9.5f %12.5f %12.5f %6.1f %6.3f %2i %6i %6.3f \n'%(tuple(map(positiveRadians,angles))+(y*1e6,x*1e6)+info[validData]) return '%9.5f %9.5f %9.5f %12.5f %12.5f %6.1f %6.3f %2i %6i %6.3f \n'%(tuple(map(positiveRadians,angles))+(y*1e6,x*1e6)+info[validData])
@ -155,10 +151,9 @@ if options.hexagonal:
for filename in filenames: for filename in filenames:
# Read the source file
if options.verbose: sys.stdout.write("\nREADING VTK FILE\n") if options.verbose: sys.stdout.write("\nREADING VTK FILE\n")
reader = vtkUnstructuredGridReader() # Read the source file
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(filename) reader.SetFileName(filename)
reader.ReadAllScalarsOn() reader.ReadAllScalarsOn()
reader.ReadAllVectorsOn() reader.ReadAllVectorsOn()
@ -166,7 +161,7 @@ for filename in filenames:
undeformedMesh = reader.GetOutput() undeformedMesh = reader.GetOutput()
# Get euler angles from cell data # Get euler angles from cell data
if options.verbose: sys.stdout.write("\nGETTING EULER ANGLES\n") if options.verbose: sys.stdout.write("\nGETTING EULER ANGLES\n")
angles = {} angles = {}
@ -177,14 +172,14 @@ for filename in filenames:
if options.verbose: sys.stdout.write(" found scalar with name %s\n"%scalarName) if options.verbose: sys.stdout.write(" found scalar with name %s\n"%scalarName)
if len(angles) < 3: # found data for all three euler angles? if len(angles) < 3: # found data for all three euler angles?
for label in options.eulerLabel: for label in options.eulerLabel:
if not label in angles.keys(): if label not in angles.keys():
parser.error('Could not find scalar data with name %s'%label) parser.error('Could not find scalar data with name %s'%label)
# Get deformed mesh # Get deformed mesh
if options.verbose: sys.stdout.write("\nDEFORM MESH\n") if options.verbose: sys.stdout.write("\nDEFORM MESH\n")
warpVector = vtkWarpVector() warpVector = vtk.vtkWarpVector()
undeformedMesh.GetPointData().SetActiveVectors(options.dispLabel) undeformedMesh.GetPointData().SetActiveVectors(options.dispLabel)
warpVector.SetInput(undeformedMesh) warpVector.SetInput(undeformedMesh)
warpVector.Update() warpVector.Update()
@ -197,29 +192,29 @@ for filename in filenames:
sys.stdout.write(" z (% .8f % .8f)\n"%(box[4],box[5])) sys.stdout.write(" z (% .8f % .8f)\n"%(box[4],box[5]))
# Get cell centers of deformed mesh (position of ips) # Get cell centers of deformed mesh (position of ips)
if options.verbose: sys.stdout.write("\nGETTING CELL CENTERS OF DEFORMED MESH\n") if options.verbose: sys.stdout.write("\nGETTING CELL CENTERS OF DEFORMED MESH\n")
cellCenter = vtkCellCenters() cellCenter = vtk.vtkCellCenters()
cellCenter.SetVertexCells(0) # do not generate vertex cells, just points cellCenter.SetVertexCells(0) # do not generate vertex cells, just points
cellCenter.SetInput(deformedMesh) cellCenter.SetInput(deformedMesh)
cellCenter.Update() cellCenter.Update()
meshIPs = cellCenter.GetOutput() meshIPs = cellCenter.GetOutput()
# Get outer surface of deformed mesh # Get outer surface of deformed mesh
if options.verbose: sys.stdout.write("\nGETTING OUTER SURFACE OF DEFORMED MESH\n") if options.verbose: sys.stdout.write("\nGETTING OUTER SURFACE OF DEFORMED MESH\n")
surfaceFilter = vtkDataSetSurfaceFilter() surfaceFilter = vtk.vtkDataSetSurfaceFilter()
surfaceFilter.SetInput(deformedMesh) surfaceFilter.SetInput(deformedMesh)
surfaceFilter.Update() surfaceFilter.Update()
surface = surfaceFilter.GetOutput() surface = surfaceFilter.GetOutput()
# Get coordinate system for ang files # Get coordinate system for ang files
# z-vector is normal to slices # z-vector is normal to slices
# x-vector corresponds to the up-direction # x-vector corresponds to the up-direction
# "R" rotates coordinates from the mesh system into the TSL system # "R" rotates coordinates from the mesh system into the TSL system
if options.verbose: sys.stdout.write("\nGETTING COORDINATE SYSTEM FOR ANG FILES\n") if options.verbose: sys.stdout.write("\nGETTING COORDINATE SYSTEM FOR ANG FILES\n")
z = np.array(options.normal,dtype='float') z = np.array(options.normal,dtype='float')
@ -235,7 +230,7 @@ for filename in filenames:
sys.stdout.write(" z (% .8f % .8f % .8f)\n"%tuple(z)) sys.stdout.write(" z (% .8f % .8f % .8f)\n"%tuple(z))
# Get bounding box in rotated system (x,y,z) # Get bounding box in rotated system (x,y,z)
if options.verbose: sys.stdout.write("\nGETTING BOUNDING BOX IN ROTATED SYSTEM\n") if options.verbose: sys.stdout.write("\nGETTING BOUNDING BOX IN ROTATED SYSTEM\n")
rotatedbox = [[np.inf,-np.inf] for i in range(3)] # bounding box in rotated TSL system rotatedbox = [[np.inf,-np.inf] for i in range(3)] # bounding box in rotated TSL system
@ -254,8 +249,8 @@ for filename in filenames:
sys.stdout.write(" z (% .8f % .8f)\n"%tuple(rotatedbox[2])) sys.stdout.write(" z (% .8f % .8f)\n"%tuple(rotatedbox[2]))
# Correct bounding box so that a multiplicity of the resolution fits into it # Correct bounding box so that a multiplicity of the resolution fits into it
# and get number of points and extent in each (rotated) axis direction # and get number of points and extent in each (rotated) axis direction
if options.verbose: sys.stdout.write("\nCORRECTING EXTENT OF BOUNDING BOX IN ROTATED SYSTEM\n") if options.verbose: sys.stdout.write("\nCORRECTING EXTENT OF BOUNDING BOX IN ROTATED SYSTEM\n")
correction = [] correction = []
@ -284,12 +279,12 @@ for filename in filenames:
sys.stdout.write(" z (% .8f % .8f)\n"%tuple(rotatedbox[2])) sys.stdout.write(" z (% .8f % .8f)\n"%tuple(rotatedbox[2]))
# Generate new regular point grid for ang files # Generate new regular point grid for ang files
# Use "polydata" object with points as single vertices # Use "polydata" object with points as single vertices
# beware of TSL convention: y direction is fastest index # beware of TSL convention: y direction is fastest index
if options.verbose: sys.stdout.write("\nGENERATING POINTS FOR POINT GRID") if options.verbose: sys.stdout.write("\nGENERATING POINTS FOR POINT GRID")
points = vtkPoints() points = vtk.vtkPoints()
for k in xrange(Npoints[2]): for k in xrange(Npoints[2]):
for j in xrange(Npoints[0]): for j in xrange(Npoints[0]):
for i in xrange(Npoints[1]): # y is fastest index for i in xrange(Npoints[1]): # y is fastest index
@ -309,9 +304,9 @@ for filename in filenames:
sys.stdout.write(" grid resolution: %.8f\n"%options.resolution) sys.stdout.write(" grid resolution: %.8f\n"%options.resolution)
if options.verbose: sys.stdout.write("\nGENERATING VERTICES FOR POINT GRID") if options.verbose: sys.stdout.write("\nGENERATING VERTICES FOR POINT GRID")
vertices = vtkCellArray() vertices = vtk.vtkCellArray()
for i in xrange(totalNpoints): for i in xrange(totalNpoints):
vertex = vtkVertex() vertex = vtk.vtkVertex()
vertex.GetPointIds().SetId(0,i) # each vertex consists of exactly one (index 0) point with ID "i" vertex.GetPointIds().SetId(0,i) # each vertex consists of exactly one (index 0) point with ID "i"
vertices.InsertNextCell(vertex) vertices.InsertNextCell(vertex)
if options.verbose: if options.verbose:
@ -319,34 +314,35 @@ for filename in filenames:
sys.stdout.flush() sys.stdout.flush()
if options.verbose: sys.stdout.write("\n\nGENERATING POINT GRID\n") if options.verbose: sys.stdout.write("\n\nGENERATING POINT GRID\n")
pointgrid = vtkPolyData() pointgrid = vtk.vtkPolyData()
pointgrid.SetPoints(points) pointgrid.SetPoints(points)
pointgrid.SetVerts(vertices) pointgrid.SetVerts(vertices)
pointgrid.Update() pointgrid.Update()
# Find out which points reside inside mesh geometry # Find out which points reside inside mesh geometry
if options.verbose: sys.stdout.write("\nIDENTIFYING POINTS INSIDE MESH GEOMETRY\n") if options.verbose: sys.stdout.write("\nIDENTIFYING POINTS INSIDE MESH GEOMETRY\n")
enclosedPoints = vtkSelectEnclosedPoints() enclosedPoints = vtk.vtkSelectEnclosedPoints()
enclosedPoints.SetSurface(surface) enclosedPoints.SetSurface(surface)
enclosedPoints.SetInput(pointgrid) enclosedPoints.SetInput(pointgrid)
enclosedPoints.Update() enclosedPoints.Update()
# Build kdtree from mesh IPs and match mesh IPs to point grid # Build kdtree from mesh IPs and match mesh IPs to point grid
if options.verbose: sys.stdout.write("\nBUILDING MAPPING OF GRID POINTS") if options.verbose: sys.stdout.write("\nBUILDING MAPPING OF GRID POINTS")
kdTree = vtkKdTree() kdTree = vtk.vtkKdTree()
kdTree.BuildLocatorFromPoints(meshIPs.GetPoints()) kdTree.BuildLocatorFromPoints(meshIPs.GetPoints())
gridToMesh = [] gridToMesh = []
ids = vtkIdList() ids = vtk.vtkIdList()
NenclosedPoints = 0 NenclosedPoints = 0
for i in range(pointgrid.GetNumberOfPoints()): for i in range(pointgrid.GetNumberOfPoints()):
gridToMesh.append([]) gridToMesh.append([])
if enclosedPoints.IsInside(i): if enclosedPoints.IsInside(i):
NenclosedPoints += 1 NenclosedPoints += 1
kdTree.FindClosestNPoints(options.interpolation,pointgrid.GetPoint(i),ids) # here one could use faster(?) "FindClosestPoint" if only first nearest neighbor required # here one could use faster(?) "FindClosestPoint" if only first nearest neighbor required
kdTree.FindClosestNPoints(options.interpolation,pointgrid.GetPoint(i),ids)
for j in range(ids.GetNumberOfIds()): for j in range(ids.GetNumberOfIds()):
gridToMesh[-1].extend([ids.GetId(j)]) gridToMesh[-1].extend([ids.GetId(j)])
if options.verbose: if options.verbose:
@ -358,7 +354,7 @@ for filename in filenames:
# ITERATE OVER SLICES AND CREATE ANG FILE # ITERATE OVER SLICES AND CREATE ANG FILE
if options.verbose: if options.verbose:
sys.stdout.write("\nWRITING OUT ANG FILES\n") sys.stdout.write("\nWRITING OUT ANG FILES\n")
@ -404,13 +400,13 @@ for filename in filenames:
angfile.write(getDataLine(interpolatedPhi,x,y,enclosedPoints.IsInside(i))) angfile.write(getDataLine(interpolatedPhi,x,y,enclosedPoints.IsInside(i)))
# Visualize slices # Visualize slices
if options.visualize: if options.visualize:
meshMapper = vtkDataSetMapper() meshMapper = vtk.vtkDataSetMapper()
meshMapper.SetInput(surface) meshMapper.SetInput(surface)
meshMapper.ScalarVisibilityOff() # do not use scalar data for coloring meshMapper.ScalarVisibilityOff() # do not use scalar data for coloring
meshActor = vtkActor() meshActor = vtk.vtkActor()
meshActor.SetMapper(meshMapper) meshActor.SetMapper(meshMapper)
meshActor.GetProperty().SetOpacity(0.2) meshActor.GetProperty().SetOpacity(0.2)
meshActor.GetProperty().SetColor(1.0,1.0,0) meshActor.GetProperty().SetColor(1.0,1.0,0)
@ -418,43 +414,43 @@ for filename in filenames:
# meshActor.GetProperty().SetEdgeColor(1,1,0.5) # meshActor.GetProperty().SetEdgeColor(1,1,0.5)
# meshActor.GetProperty().EdgeVisibilityOn() # meshActor.GetProperty().EdgeVisibilityOn()
boxpoints = vtkPoints() boxpoints = vtk.vtkPoints()
for n in range(8): for n in range(8):
P = [rotatedbox[0][(n/1)%2], P = [rotatedbox[0][(n/1)%2],
rotatedbox[1][(n/2)%2], rotatedbox[1][(n/2)%2],
rotatedbox[2][(n/4)%2]] rotatedbox[2][(n/4)%2]]
boxpoints.InsertNextPoint(list(np.dot(R.T,np.array(P)))) boxpoints.InsertNextPoint(list(np.dot(R.T,np.array(P))))
box = vtkHexahedron() box = vtk.vtkHexahedron()
for n,i in enumerate([0,1,3,2,4,5,7,6]): for n,i in enumerate([0,1,3,2,4,5,7,6]):
box.GetPointIds().SetId(n,i) box.GetPointIds().SetId(n,i)
boxgrid = vtkUnstructuredGrid() boxgrid = vtk.vtkUnstructuredGrid()
boxgrid.SetPoints(boxpoints) boxgrid.SetPoints(boxpoints)
boxgrid.InsertNextCell(box.GetCellType(), box.GetPointIds()) boxgrid.InsertNextCell(box.GetCellType(), box.GetPointIds())
boxsurfaceFilter = vtkDataSetSurfaceFilter() boxsurfaceFilter = vtk.vtkDataSetSurfaceFilter()
boxsurfaceFilter.SetInput(boxgrid) boxsurfaceFilter.SetInput(boxgrid)
boxsurfaceFilter.Update() boxsurfaceFilter.Update()
boxsurface = boxsurfaceFilter.GetOutput() boxsurface = boxsurfaceFilter.GetOutput()
boxMapper = vtkDataSetMapper() boxMapper = vtk.vtkDataSetMapper()
boxMapper.SetInput(boxsurface) boxMapper.SetInput(boxsurface)
boxActor = vtkActor() boxActor = vtk.vtkActor()
boxActor.SetMapper(boxMapper) boxActor.SetMapper(boxMapper)
boxActor.GetProperty().SetLineWidth(2.0) boxActor.GetProperty().SetLineWidth(2.0)
boxActor.GetProperty().SetRepresentationToWireframe() boxActor.GetProperty().SetRepresentationToWireframe()
gridMapper = vtkDataSetMapper() gridMapper = vtk.vtkDataSetMapper()
gridMapper.SetInput(pointgrid) gridMapper.SetInput(pointgrid)
gridActor = vtkActor() gridActor = vtk.vtkActor()
gridActor.SetMapper(gridMapper) gridActor.SetMapper(gridMapper)
gridActor.GetProperty().SetColor(0,0,0) gridActor.GetProperty().SetColor(0,0,0)
gridActor.GetProperty().SetPointSize(3) gridActor.GetProperty().SetPointSize(3)
renderer = vtkRenderer() renderer = vtk.vtkRenderer()
renderWindow = vtkRenderWindow() renderWindow = vtk.vtkRenderWindow()
renderWindow.FullScreenOn() renderWindow.FullScreenOn()
renderWindow.AddRenderer(renderer) renderWindow.AddRenderer(renderer)
renderWindowInteractor = vtkRenderWindowInteractor() renderWindowInteractor = vtk.vtkRenderWindowInteractor()
renderWindowInteractor.SetRenderWindow(renderWindow) renderWindowInteractor.SetRenderWindow(renderWindow)
renderer.AddActor(meshActor) renderer.AddActor(meshActor)
renderer.AddActor(boxActor) renderer.AddActor(boxActor)
@ -462,6 +458,6 @@ for filename in filenames:
renderer.SetBackground(1,1,1) renderer.SetBackground(1,1,1)
renderWindow.Render() renderWindow.Render()
renderWindowInteractor.SetInteractorStyle(vtkInteractorStyleTrackballCamera()) renderWindowInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
renderWindowInteractor.Start() renderWindowInteractor.Start()

1
processing/post/vtk_addData.py
View File

@ -10,7 +10,6 @@ scriptID = ' '.join([scriptName,damask.version])
# ----------------------------- # -----------------------------
def findTag(filename,tag): def findTag(filename,tag):
# -----------------------------
with open(filename,'r') as myfile: with open(filename,'r') as myfile:
mypattern = re.compile(str(tag)) mypattern = re.compile(str(tag))

20
processing/post/vtk_scalars2vectors.py
View File

@ -4,7 +4,7 @@
import os,sys,shutil import os,sys,shutil
import damask import damask
from optparse import OptionParser from optparse import OptionParser
from vtk import * import vtk
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])
@ -44,18 +44,16 @@ for filename in filenames:
for filename in filenames: for filename in filenames:
# Read the source file
sys.stdout.write('read file "%s" ...'%filename) sys.stdout.write('read file "%s" ...'%filename)
sys.stdout.flush() sys.stdout.flush()
suffix = os.path.splitext(filename)[1] suffix = os.path.splitext(filename)[1]
if suffix == '.vtk': if suffix == '.vtk':
reader = vtkUnstructuredGridReader() reader = vtk.vtkUnstructuredGridReader()
reader.ReadAllScalarsOn() reader.ReadAllScalarsOn()
reader.ReadAllVectorsOn() reader.ReadAllVectorsOn()
reader.ReadAllTensorsOn() reader.ReadAllTensorsOn()
elif suffix == '.vtu': elif suffix == '.vtu':
reader = vtkXMLUnstructuredGridReader() reader = vtk.vtkXMLUnstructuredGridReader()
else: else:
parser.error('filetype "%s" not supported'%suffix) parser.error('filetype "%s" not supported'%suffix)
reader.SetFileName(filename) reader.SetFileName(filename)
@ -65,7 +63,7 @@ for filename in filenames:
sys.stdout.flush() sys.stdout.flush()
# Read the scalar data # Read the scalar data
scalarData = {} scalarData = {}
scalarsToBeRemoved = [] scalarsToBeRemoved = []
@ -83,19 +81,18 @@ for filename in filenames:
scalarsToBeRemoved.append(scalarName) scalarsToBeRemoved.append(scalarName)
for scalarName in scalarsToBeRemoved: for scalarName in scalarsToBeRemoved:
uGrid.GetCellData().RemoveArray(scalarName) uGrid.GetCellData().RemoveArray(scalarName)
# uGrid.UpdateData()
sys.stdout.write('\rread scalar data done\n') sys.stdout.write('\rread scalar data done\n')
sys.stdout.flush() sys.stdout.flush()
# Convert the scalar data to vector data # Convert the scalar data to vector data
NscalarData = len(scalarData) NscalarData = len(scalarData)
for n,label in enumerate(scalarData): for n,label in enumerate(scalarData):
sys.stdout.write("\rconvert to vector data %d%%" %(100*n/NscalarData)) sys.stdout.write("\rconvert to vector data %d%%" %(100*n/NscalarData))
sys.stdout.flush() sys.stdout.flush()
Nvalues = scalarData[label][0].GetNumberOfTuples() Nvalues = scalarData[label][0].GetNumberOfTuples()
vectorData = vtkDoubleArray() vectorData = vtk.vtkDoubleArray()
vectorData.SetName(label) vectorData.SetName(label)
vectorData.SetNumberOfComponents(3) # set this before NumberOfTuples !!! vectorData.SetNumberOfComponents(3) # set this before NumberOfTuples !!!
vectorData.SetNumberOfTuples(Nvalues) vectorData.SetNumberOfTuples(Nvalues)
@ -103,16 +100,15 @@ for filename in filenames:
for j in range(3): for j in range(3):
vectorData.SetComponent(i,j,scalarData[label][j].GetValue(i)) vectorData.SetComponent(i,j,scalarData[label][j].GetValue(i))
uGrid.GetCellData().AddArray(vectorData) uGrid.GetCellData().AddArray(vectorData)
# uGrid.GetCellData().SetActiveVectors(label)
sys.stdout.write('\rconvert to vector data done\n') sys.stdout.write('\rconvert to vector data done\n')
# Write to new vtk file # Write to new vtk file
outfilename = os.path.splitext(filename)[0]+'.vtu' outfilename = os.path.splitext(filename)[0]+'.vtu'
sys.stdout.write('write to file "%s" ...'%outfilename) sys.stdout.write('write to file "%s" ...'%outfilename)
sys.stdout.flush() sys.stdout.flush()
writer = vtkXMLUnstructuredGridWriter() writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName(outfilename+'_tmp') writer.SetFileName(outfilename+'_tmp')
writer.SetDataModeToAscii() writer.SetDataModeToAscii()
writer.SetInput(uGrid) writer.SetInput(uGrid)