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

WIP: making compatible with python3/vtk9

This commit is contained in:
Martin Diehl 2020-08-08 18:19:04 +02:00
parent 0ad189ea9d
commit 0878302961
1 changed files with 56 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

130
processing/post/vtk2ang.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7 #!/usr/bin/env python3
import os import os
import sys import sys
@ -11,13 +11,13 @@ import vtk
import damask 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', \
@ -53,7 +53,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 from 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
@ -67,7 +67,7 @@ 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 Convention in ang file: y coordinate comes first and is fastest index
positions in micrometer 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)}
@ -75,18 +75,16 @@ def getDataLine(angles,x,y,validData=True):
%(tuple(map(positiveRadians,angles))+(y*1e6,x*1e6)+info[validData]) %(tuple(map(positiveRadians,angles))+(y*1e6,x*1e6)+info[validData])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN FUNCTION STARTS HERE # MAIN FUNCTION STARTS HERE
# -------------------------------------------------------------------- # --------------------------------------------------------------------
parser = OptionParser(usage='%prog options [file[s]]', description = """ parser = OptionParser(usage='%prog options [file[s]]', description = """
Builds a ang files from a vtk file. Builds a ang files from a vtk file.
""", version = scriptID) """, version = scriptID)
parser.add_option('--disp','--displacement',dest='dispLabel', parser.add_option('--disp','--displacement',dest='dispLabel',
metavar ='string', metavar ='string',
help='label of displacements [%default]') help='label of displacements [%default]')
parser.add_option('--euler', dest='eulerLabel', nargs=3, parser.add_option('--euler', dest='eulerLabel', nargs=3,
@ -110,11 +108,6 @@ parser.add_option('-s','--scale', dest='scale', type='float',
parser.add_option('-r','--resolution', dest='resolution', type='float', parser.add_option('-r','--resolution', dest='resolution', type='float',
metavar ='float', metavar ='float',
help='scaling factor for resolution [%default]') help='scaling factor for resolution [%default]')
parser.add_option('--hex','--hexagonal', dest='hexagonal', action='store_true',
help='use in plane hexagonal grid')
parser.add_option('--interpolation', dest='interpolation', type='int',
metavar='float',
help='number of points for linear interpolation [%default]')
parser.add_option('--verbose', dest='verbose', action='store_true', parser.add_option('--verbose', dest='verbose', action='store_true',
help='verbose mode') help='verbose mode')
parser.add_option('--visualize', dest='visualize', action='store_true', parser.add_option('--visualize', dest='visualize', action='store_true',
@ -122,7 +115,6 @@ parser.add_option('--visualize', dest='visualize', action='store_true
parser.set_defaults(dispLabel = 'displacement') parser.set_defaults(dispLabel = 'displacement')
parser.set_defaults(eulerLabel = ['1_1_eulerangles','1_2_eulerangles','1_3_eulerangles']) parser.set_defaults(eulerLabel = ['1_1_eulerangles','1_2_eulerangles','1_3_eulerangles'])
parser.set_defaults(hexagonal = False)
parser.set_defaults(normal = [0.0,0.0,-1.0]) parser.set_defaults(normal = [0.0,0.0,-1.0])
parser.set_defaults(up = [0.0,1.0,0.0]) parser.set_defaults(up = [0.0,1.0,0.0])
parser.set_defaults(Nslices = 1) parser.set_defaults(Nslices = 1)
@ -130,7 +122,6 @@ parser.set_defaults(distance = 0.0)
parser.set_defaults(scale = 1.0) parser.set_defaults(scale = 1.0)
parser.set_defaults(resolution = 1.0) parser.set_defaults(resolution = 1.0)
parser.set_defaults(dispScaling = 1.0) parser.set_defaults(dispScaling = 1.0)
parser.set_defaults(interpolation = 1)
parser.set_defaults(verbose = False) parser.set_defaults(verbose = False)
parser.set_defaults(visualize = False) parser.set_defaults(visualize = False)
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -153,35 +144,26 @@ if np.dot(np.array(options.normal),np.array(options.up)) > 1e-3:
parser.error('normal vector and up vector have to be orthogonal') parser.error('normal vector and up vector have to be orthogonal')
# check for options that are not yet implemented
if options.interpolation > 1:
parser.error('interpolation not yet supported')
if options.hexagonal:
parser.error('hexagonal grid not yet supported')
#--- ITERATE OVER FILES AND PROCESS THEM #--- ITERATE OVER FILES AND PROCESS THEM
for filename in filenames: for filename in filenames:
if options.verbose: sys.stdout.write("\nREADING VTK FILE\n") if options.verbose: sys.stdout.write("\nREADING VTK FILE\n")
# Read the source file # Read the source file
reader = vtk.vtkUnstructuredGridReader() reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(filename) reader.SetFileName(filename)
reader.ReadAllScalarsOn() reader.ReadAllScalarsOn()
reader.ReadAllVectorsOn() reader.ReadAllVectorsOn()
reader.Update() reader.Update()
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 = {}
for i in range(reader.GetNumberOfScalarsInFile()): for i in range(undeformedMesh.GetPointData().GetNumberOfArrays()):
scalarName = reader.GetScalarsNameInFile(i) scalarName = undeformedMesh.GetPointData().GetArrayName(i)
if scalarName in options.eulerLabel: if scalarName in options.eulerLabel:
angles[scalarName] = undeformedMesh.GetCellData().GetScalars(scalarName) angles[scalarName] = undeformedMesh.GetCellData().GetScalars(scalarName)
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)
@ -189,14 +171,14 @@ for filename in filenames:
for label in options.eulerLabel: for label in options.eulerLabel:
if label not 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 = vtk.vtkWarpVector() warpVector = vtk.vtkWarpVector()
undeformedMesh.GetPointData().SetActiveVectors(options.dispLabel) undeformedMesh.GetPointData().SetActiveVectors(options.dispLabel)
warpVector.SetInput(undeformedMesh) warpVector.SetInputData(undeformedMesh)
warpVector.Update() warpVector.Update()
deformedMesh = warpVector.GetOutput() deformedMesh = warpVector.GetOutput()
box = deformedMesh.GetBounds() # bounding box in mesh system box = deformedMesh.GetBounds() # bounding box in mesh system
@ -208,24 +190,24 @@ for filename in filenames:
# 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 = vtk.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.SetInputData(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 = vtk.vtkDataSetSurfaceFilter() surfaceFilter = vtk.vtkDataSetSurfaceFilter()
surfaceFilter.SetInput(deformedMesh) surfaceFilter.SetInputData(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
@ -249,10 +231,10 @@ for filename in filenames:
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
for n in range(8): # loop over eight vertices of mesh bounding box for n in range(8): # loop over eight vertices of mesh bounding box
vert = np.array([box[0+(n/1)%2], vert = np.array([box[0+(n//1)%2],
box[2+(n/2)%2], box[2+(n//2)%2],
box[4+(n/4)%2]]) # vertex in mesh system box[4+(n//4)%2]]) # vertex in mesh system
rotatedvert = np.dot(R,vert) # vertex in rotated system rotatedvert = np.dot(R,vert) # vertex in rotated system
for i in range(3): for i in range(3):
rotatedbox[i][0] = min(rotatedbox[i][0],rotatedvert[i]) rotatedbox[i][0] = min(rotatedbox[i][0],rotatedvert[i])
@ -274,7 +256,7 @@ for filename in filenames:
for i in range(2): for i in range(2):
Npoints.extend([int(math.ceil(extent[i] / options.resolution))]) Npoints.extend([int(math.ceil(extent[i] / options.resolution))])
correction.extend([float(Npoints[i]) * options.resolution - extent[i]]) correction.extend([float(Npoints[i]) * options.resolution - extent[i]])
if options.distance > 0.0: if options.distance > 0.0:
Npoints.extend([int(math.ceil(extent[2] / options.distance))]) Npoints.extend([int(math.ceil(extent[2] / options.distance))])
correction.extend([float(Npoints[2]) * options.distance - extent[2]]) correction.extend([float(Npoints[2]) * options.distance - extent[2]])
else: else:
@ -282,8 +264,8 @@ for filename in filenames:
correction.extend([0.0]) correction.extend([0.0])
options.distance = extent[2] / float(options.Nslices) options.distance = extent[2] / float(options.Nslices)
for i in range(3): for i in range(3):
rotatedbox[i][0] = rotatedbox[i][0] - 0.5 * correction[i] rotatedbox[i][0] -= 0.5 * correction[i]
rotatedbox[i][1] = rotatedbox[i][1] + 0.5 * correction[i] rotatedbox[i][1] += 0.5 * correction[i]
extent[i] = rotatedbox[i][1] - rotatedbox[i][0] extent[i] = rotatedbox[i][1] - rotatedbox[i][0]
NpointsPerSlice = Npoints[0] * Npoints[1] NpointsPerSlice = Npoints[0] * Npoints[1]
totalNpoints = NpointsPerSlice * Npoints[2] totalNpoints = NpointsPerSlice * Npoints[2]
@ -301,20 +283,20 @@ for filename in filenames:
if options.verbose: sys.stdout.write("\nGENERATING POINTS FOR POINT GRID") if options.verbose: sys.stdout.write("\nGENERATING POINTS FOR POINT GRID")
points = vtk.vtkPoints() points = vtk.vtkPoints()
for k in range(Npoints[2]): for k in range(Npoints[2]):
for j in range(Npoints[0]): for j in range(Npoints[0]):
for i in range(Npoints[1]): # y is fastest index for i in range(Npoints[1]): # y is fastest index
rotatedpoint = np.array([rotatedbox[0][0] + (float(j) + 0.5) * options.resolution, rotatedpoint = np.array([rotatedbox[0][0] + (float(j) + 0.5) * options.resolution,
rotatedbox[1][0] + (float(i) + 0.5) * options.resolution, rotatedbox[1][0] + (float(i) + 0.5) * options.resolution,
rotatedbox[2][0] + (float(k) + 0.5) * options.distance ]) # point in rotated system rotatedbox[2][0] + (float(k) + 0.5) * options.distance ]) # point in rotated system
point = np.dot(R.T,rotatedpoint) # point in mesh system point = np.dot(R.T,rotatedpoint) # point in mesh system
points.InsertNextPoint(list(point)) points.InsertNextPoint(list(point))
if options.verbose: if options.verbose:
sys.stdout.write("\rGENERATING POINTS FOR POINT GRID %d%%" %(100*(Npoints[1]*(k*Npoints[0]+j)+i+1)/totalNpoints)) sys.stdout.write("\rGENERATING POINTS FOR POINT GRID %d%%" %(100*(Npoints[1]*(k*Npoints[0]+j)+i+1)/totalNpoints))
sys.stdout.flush() sys.stdout.flush()
if options.verbose: if options.verbose:
sys.stdout.write("\n number of slices: %i\n"%Npoints[2]) sys.stdout.write("\n number of slices: %i\n"%Npoints[2])
sys.stdout.write(" slice spacing: %.8f\n"%options.distance) sys.stdout.write(" slice spacing: %.8f\n"%options.distance)
if Npoints[2] > 1: if Npoints[2] > 1:
sys.stdout.write(" number of points per slice: %i = %i rows * %i points in row\n"%(NpointsPerSlice,Npoints[0],Npoints[1])) sys.stdout.write(" number of points per slice: %i = %i rows * %i points in row\n"%(NpointsPerSlice,Npoints[0],Npoints[1]))
sys.stdout.write(" grid resolution: %.8f\n"%options.resolution) sys.stdout.write(" grid resolution: %.8f\n"%options.resolution)
@ -324,7 +306,7 @@ for filename in filenames:
vertex = vtk.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:
sys.stdout.write("\rGENERATING VERTICES FOR POINT GRID %d%%" %(100*(i+1)/totalNpoints)) sys.stdout.write("\rGENERATING VERTICES FOR POINT GRID %d%%" %(100*(i+1)/totalNpoints))
sys.stdout.flush() sys.stdout.flush()
@ -357,34 +339,34 @@ for filename in filenames:
if enclosedPoints.IsInside(i): if enclosedPoints.IsInside(i):
NenclosedPoints += 1 NenclosedPoints += 1
# 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) kdTree.FindClosestNPoints(1,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:
sys.stdout.write("\rBUILDING MAPPING OF GRID POINTS %d%%" %(100*(i+1)/totalNpoints)) sys.stdout.write("\rBUILDING MAPPING OF GRID POINTS %d%%" %(100*(i+1)/totalNpoints))
sys.stdout.flush() sys.stdout.flush()
if options.verbose: if options.verbose:
sys.stdout.write("\n Number of points inside mesh geometry %i\n"%NenclosedPoints) sys.stdout.write("\n Number of points inside mesh geometry %i\n"%NenclosedPoints)
sys.stdout.write(" Number of points outside mesh geometry %i\n"%(totalNpoints - NenclosedPoints)) sys.stdout.write(" Number of points outside mesh geometry %i\n"%(totalNpoints - NenclosedPoints))
# 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")
sys.stdout.write(" scaling all length with %f\n"%options.scale) sys.stdout.write(" scaling all length with %f\n"%options.scale)
x0,y0,z0 = np.dot(R,pointgrid.GetPoint(0)) # first point on slice defines origin x0,y0,z0 = np.dot(R,pointgrid.GetPoint(0)) # first point on slice defines origin
for sliceN in range(Npoints[2]): for sliceN in range(Npoints[2]):
# Open file and write header # Open file and write header
angfilename = eval('"'+eval("'%%s_slice%%0%ii.ang'%(math.log10(Npoints[2])+1)")+'"%(os.path.splitext(filename)[0],sliceN+1)') angfilename = eval('"'+eval("'%%s_slice%%0%ii.ang'%(math.log10(Npoints[2])+1)")+'"%(os.path.splitext(filename)[0],sliceN+1)')
with open(angfilename,'w') as angfile: with open(angfilename,'w') as angfile:
if options.verbose: sys.stdout.write(" %s\n"%angfilename) if options.verbose: sys.stdout.write(" %s\n"%angfilename)
angfile.write(getHeader(filename,Npoints[1],Npoints[0],options.resolution*options.scale)) angfile.write(getHeader(filename,Npoints[1],Npoints[0],options.resolution*options.scale))
for i in range(sliceN*NpointsPerSlice,(sliceN+1)*NpointsPerSlice): # Iterate over points on slice for i in range(sliceN*NpointsPerSlice,(sliceN+1)*NpointsPerSlice): # Iterate over points on slice
# Get euler angles of closest IDs # Get euler angles of closest IDs
@ -397,26 +379,26 @@ for filename in filenames:
phi[-1].extend([angles[options.eulerLabel[k]].GetValue(IP)]) phi[-1].extend([angles[options.eulerLabel[k]].GetValue(IP)])
else: else:
phi = [[720,720,720]] # fake angles phi = [[720,720,720]] # fake angles
# Interpolate Euler angle # Interpolate Euler angle
# NOT YET IMPLEMENTED, simply take the nearest neighbors values # NOT YET IMPLEMENTED, simply take the nearest neighbors values
interpolatedPhi = phi[0] interpolatedPhi = phi[0]
# write data to ang file # write data to ang file
x,y,z = np.dot(R,pointgrid.GetPoint(i)) # point in rotated TSL system x,y,z = np.dot(R,pointgrid.GetPoint(i)) # point in rotated TSL system
x -= x0 # first point on slice defines origin x -= x0 # first point on slice defines origin
y -= y0 # first point on slice defines origin y -= y0 # first point on slice defines origin
x *= options.scale x *= options.scale
y *= options.scale y *= options.scale
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 = vtk.vtkDataSetMapper() meshMapper = vtk.vtkDataSetMapper()
meshMapper.SetInput(surface) meshMapper.SetInput(surface)
@ -426,11 +408,11 @@ for filename in filenames:
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)
meshActor.GetProperty().BackfaceCullingOn() meshActor.GetProperty().BackfaceCullingOn()
boxpoints = vtk.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 = vtk.vtkHexahedron() box = vtk.vtkHexahedron()
@ -469,7 +451,7 @@ for filename in filenames:
renderer.AddActor(boxActor) renderer.AddActor(boxActor)
renderer.AddActor(gridActor) renderer.AddActor(gridActor)
renderer.SetBackground(1,1,1) renderer.SetBackground(1,1,1)
renderWindow.Render() renderWindow.Render()
renderWindowInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera()) renderWindowInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
renderWindowInteractor.Start() renderWindowInteractor.Start()