Skip to content
Snippets Groups Projects
Commit ecc5d3f5 authored by David Seus's avatar David Seus
Browse files

remove Jim's python plotscript due to copyright concerns

parent 439d9e3a
No related branches found
No related tags found
No related merge requests found
https://lorensen.github.io/VTKExamples/site/Python/
#!/usr/bin/python
# Python script to plot .vtk and .vtu data
# and adding streamlines.
import os
from os.path import basename
import sys
from optparse import OptionParser
import subprocess
import re
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.tri as mtri
from matplotlib import cm
import matplotlib.patches
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.colors import LinearSegmentedColormap
import fileinput
from glob import glob
from scipy.interpolate import griddata
import numpy as np
from vtk import *
from vtk.util.numpy_support import vtk_to_numpy
from vtk import vtkUnstructuredGridReader
from vtk import vtkXMLUnstructuredGridReader
def cm2inch(*tupl):
# Convert centimeters to inch
inch = 2.54
if isinstance(tupl[0], tuple):
return tuple(i / inch for i in tupl[0])
else:
return tuple(i / inch for i in tupl)
quiver_pos_x = None
quiver_pos_y = None
#========================================================================
def plot_vtk_file(filename, interface_filename=None, scalar_variable=None, datarange = None, start_points = None, cmap_binedges = None, plot_velocity=False, output_folder='default', cmap=plt.cm.Spectral_r):
quivercolor = 'white'
interfacecolor = 'black'
# cmap = plt.cm.Spectral_r
# # cmap = plt.cm.cividis
# cmap = plt.cm.RdBu_r
# # cmap = plt.cm.Blues
# Get the filetype
extension = os.path.splitext(filename)[1]
if (extension == '.vtu'):
reader = vtkXMLUnstructuredGridReader()
elif extension == '.vtk':
reader = vtkUnstructuredGridReader()
reader.SetFileName(filename)
# reader.ReadAllVectorsOn()
# reader.ReadAllScalarsOn()
reader.Update()
if interface_filename is not None:
# Get the filetype
interface_extension = os.path.splitext(interface_filename)[1]
if (interface_extension == '.vtu'):
interface_reader = vtkXMLUnstructuredGridReader()
elif interface_extension == '.vtk':
interface_reader = vtkUnstructuredGridReader()
interface_reader.SetFileName(interface_filename)
interface_reader.Update()
# Get the coordinates of nodes in the mesh
nodes_vtk_array = reader.GetOutput().GetPoints().GetData()
print('--- Read cell values...')
n_cell_arrays = reader.GetOutput().GetCellData().GetNumberOfArrays()
n_point_arrays = reader.GetOutput().GetPointData().GetNumberOfArrays()
scalar_idx = -1
for i in range(n_cell_arrays):
var_name = reader.GetOutput().GetCellData().GetArrayName(i)
if (var_name == scalar_variable):
scalar_idx = i
# print("scalar_idx = {}".format(scalar_idx))
if scalar_idx >= 0:
scalar_vtk_array = reader.GetOutput().GetCellData().GetArray(scalar_idx)
# time_field = reader.GetOutput().GetFieldData().GetArray(0)
# simulation_time = vtk_to_numpy(time_field)[0]
#####################################
nodes_numpy_array = vtk_to_numpy(nodes_vtk_array)
x, y, z = nodes_numpy_array[:, 0], nodes_numpy_array[
:, 1], nodes_numpy_array[:, 2]
if scalar_idx >= 0:
scalar_numpy_array = vtk_to_numpy(scalar_vtk_array)
if (extension == '.vtk'):
scalar_numpy_array = scalar_numpy_array[:,scalar_idx]
try:
if (np.shape(scalar_numpy_array)[1] > 1):
scalar_numpy_array = np.linalg.norm(scalar_numpy_array, axis=1)
except Exception as e:
print(e)
# raise
tri_vtk_array = reader.GetOutput().GetCells().GetData()
tri = vtk_to_numpy(tri_vtk_array)
triangles = tri
ntri = np.shape(triangles)
ntri = (ntri[0]) / (3 + 1)
ntri = int(ntri)
ia = np.zeros(ntri, dtype=int)
ib = np.zeros(ntri, dtype=int)
ic = np.zeros(ntri, dtype=int)
cell_centers_x = np.zeros(ntri)
cell_centers_y = np.zeros(ntri)
# triangles are formatted like (nPoints, id1, id2, id3, ...., idnPoints)
for i in range(0, ntri):
ia[i] = int(triangles[4 * i + 1])
ib[i] = int(triangles[4 * i + 2])
ic[i] = int(triangles[4 * i + 3])
cell_centers_x[i] = (x[ia[i]] + x[ib[i]] + x[ic[i]]) / 3.0
cell_centers_y[i] = (y[ia[i]] + y[ib[i]] + y[ic[i]]) / 3.0
triangles = np.vstack((ia, ib, ic))
triangles = triangles.T
#triangles = np.reshape(triangles, (-1, M))
# print(triangles)
# print(np.shape(triangles))
############################################
triangulation = mtri.Triangulation(x=x, y=y, triangles=triangles)
# ############################################
#
# print('--- Map to grid...')
#
# # make velocity mesh grid:
# if(plot_velocity):
# resolution = 200
# else:
# resolution = 100
# xi, yi = np.meshgrid(np.linspace(x.min(), x.max(), resolution),
# np.linspace(y.min(), y.max(), resolution))
#
# # 'cubic', 'nearest', 'linear'
# if (plot_velocity):
# ui = griddata((cell_centers_x, cell_centers_y), u, (xi, yi), method='cubic')
# vi = griddata((cell_centers_x, cell_centers_y), v, (xi, yi), method='cubic')
#
# print('--- Mask values...')
# print('--- (Note: Please adapt this manually to your application)')
# # Mask Values:
# mask = grid_mask(xi, yi)
#
# if (plot_velocity):
# vi[mask] = np.nan
# ui[mask] = np.nan
#
# vel_magnitude = np.sqrt(ui * ui + vi * vi)
#
# ############################################
print('--- Plot...')
matplotlib.style.use('ggplot')
# matplotlib.style.use('default')
#matplotlib.style.use(['seaborn-deep', 'seaborn-talk'])
# matplotlib.style.use('seaborn-darkgrid')
# matplotlib.style.use('seaborn-whitegrid')
matplotlib.rc('font', family='sans-serif')
matplotlib.rc('font', serif='Biolinum')
matplotlib.rc('text', usetex='false')
# plt.figure(figsize=(12,12))
#plt.figure(figsize=cm2inch(30, 30))
# plt.figure(figsize=cm2inch(28, 28))
plt.figure(figsize=cm2inch(15, 15))
ax = plt.gca()
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())
# plt.gca().set_aspect('equal')
ax.set_aspect('equal')
# cmap_lin = plt.cm.Spectral
# cmap = plt.cm.Spectral_r
# cmap_lin = plt.cm.RdBu
# cmap_lin = plt.get_cmap('plasma')
# cmap_lin = plt.get_cmap('viridis')
# cmap_lin = matplotlib.colors.ListedColormap(cmap_lin.colors[::-1])
if cmap_binedges is not None:
norm = matplotlib.colors.BoundaryNorm(boundaries=cmap_binedges, ncolors=256)
else:
norm = None
#plt.tripcolor(x, y, triangles, facecolors=scalar_numpy_array, edgecolors='k', cmap=plt.cm.Spectral)
# im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.cm.Spectral_r)
# im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.cm.Spectral)
if scalar_idx >= 0:
im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=cmap, norm=norm)
# im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.get_cmap('viridis'))
#im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.cm.viridis)
#plt.tripcolor(x, y, triangles, facecolors=scalar_numpy_array, cmap=plt.cm.Spectral, shading='gouraud')
#im = plt.triplot(x, y, triangles, color='grey', alpha=0.2, linewidth=0.25)
ax.triplot(triangulation, color='black', alpha=0.15, linewidth=0.2)
# plt.title('t = {:.2e}'.format(simulation_time))
plt.xlabel('x')
plt.ylabel('y')
if datarange is not None:
if (len(datarange) == 2):
plt.clim(datarange[0], datarange[1])
if interface_filename is not None:
nodes_vtk_array = interface_reader.GetOutput().GetPoints().GetData()
nodes_numpy_array = vtk_to_numpy(nodes_vtk_array)
xin, yin, zin = nodes_numpy_array[:, 0], nodes_numpy_array[:, 1], nodes_numpy_array[:, 2]
int_vtk_array = interface_reader.GetOutput().GetCells().GetData()
interface_array = vtk_to_numpy(int_vtk_array)
ninterface = np.shape(interface_array)
ninterface = (ninterface[0]) / (2 + 1)
ninterface = int(ninterface)
# cells are formatted like (nPoints, id1, id2, id3, ...., idnPoints)
for i in range(0, ninterface):
ia = int(interface_array[3 * i + 1])
ib = int(interface_array[3 * i + 2])
ax.plot([xin[ia], xin[ib]], [yin[ia], yin[ib]], alpha=1, lw=0.5, color=interfacecolor, solid_capstyle='round')
# interface_centers_x[i] = (xin[ia[i]] + xin[ib[i]]) / 2.0
# interface_centers_y[i] = (yin[ia[i]] + yin[ib[i]]) / 2.0
# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.1 inch.
if scalar_idx >= 0:
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.3)
cbar = plt.colorbar(im, cax=cax, spacing='proportional')
# cbar = plt.colorbar(im, cax=cax)
# Name the colorbar automatically:
scalar_label = reader.GetOutput().GetCellData().GetArrayName(scalar_idx)
cbar.set_label(scalar_label, labelpad=5)
# cbar.set_label(scalar_variable, labelpad=5)
if (plot_velocity):
# plot_streams = True
plot_streams = False
# plot_quiver = False
plot_quiver = True
plot_contour = False
for i in range(n_cell_arrays):
var_name = reader.GetOutput().GetCellData().GetArrayName(i)
if (var_name == 'velocity'):
vector_idx = i
# print("vector_idx = {}".format(vector_idx))
vector_vtk_array = reader.GetOutput().GetCellData().GetArray(vector_idx)
v_numpy_array = vtk_to_numpy(vector_vtk_array)
u = v_numpy_array[:,0]
v = v_numpy_array[:,1]
if(plot_streams):
print('--- add streamlines...')
#ax.streamplot(xi, yi, ui, vi, color=vel_magnitude, linewidth=1, cmap=plt.cm.YlGnBl_r)
lw = vel_magnitude
lw[np.abs(lw) < 1e-9] = 0
lw = lw / np.nanmax(lw)
if(np.nanmax(lw) != np.nanmin(lw)):
lw = 2 * lw
# cont = ax.contour(xi, yi, vel_magnitude, 3, cmap=plt.get_cmap('plasma'), linewidths=0.5, alpha=0.33)
if (start_points is None):
stream_density = 0.66
# stream_density = 0.5
stream_objects = ax.streamplot(
xi, yi, ui, vi, color=vel_magnitude, linewidth=lw, arrowsize=3.0, cmap=plt.get_cmap('plasma'), density=stream_density)
# xi, yi, ui, vi, color=vel_magnitude, linewidth=1.2, arrowsize=1.33, cmap=plt.get_cmap('plasma'), density=0.75)
else:
xxi = np.linspace(x.min(), x.max(), resolution)
yyi = np.linspace(y.min(), y.max(), resolution)
stream_objects = ax.streamplot(
xxi, yyi, ui, vi, color=vel_magnitude, linewidth=lw, arrowsize=3.0, cmap=plt.get_cmap('plasma'), start_points=start_points.T, density=200)
# xi, yi, ui, vi, color=vel_magnitude, linewidth=1.2, arrowsize=1.33, cmap=plt.get_cmap('plasma'), density=0.75)
cax_streams = divider.append_axes("right", size="3%", pad=0.7)
cbar_streams = plt.colorbar(stream_objects.lines, cax=cax_streams)
#cbar_streams.ax.set_ylabel('Velocity', rotation=270)
cbar_streams.set_label('velocity', labelpad=5)
stream_alpha = 0.8
# stream_alpha = 0.5
stream_objects.lines.set_alpha(stream_alpha)
stream_objects.arrows.set_alpha(stream_alpha)
# for x in plt.gca().get_children():
# if type(x) == matplotlib.patches.FancyArrowPatch:
# x.set_alpha(stream_alpha) # or x.set_visible(False)
# start_points = [start_points, 0.5 * (x.posA + x.posB)]
# print(start_points)
if(plot_contour):
n_contourlines = 3
cont = ax.contour(xi, yi, vel_magnitude, n_contourlines, cmap=plt.get_cmap('plasma'), linewidths=0.5, alpha=0.33)
cax_streams = divider.append_axes("right", size="3%", pad=0.7)
cbar_streams = plt.colorbar(cont, cax=cax_streams)
#cbar_streams.ax.set_ylabel('Velocity', rotation=270)
cbar_streams.set_label('velocity', labelpad=5)
if(plot_quiver):
print('Plot quiver...')
global quiver_pos_x
global quiver_pos_y
if quiver_pos_x is None:
# skip=slice(None,None,25)
skip=slice(None,None,75)
quiver_pos_x = cell_centers_x[skip]
quiver_pos_y = cell_centers_y[skip]
finder = triangulation.get_trifinder()
tri_indices = finder(quiver_pos_x, quiver_pos_y)
mu = u[tri_indices].copy()
mv = v[tri_indices].copy()
mask = np.sqrt(mu**2 + mv**2) < 1e-12
mx = np.ma.masked_where(mask, quiver_pos_x)
my = np.ma.masked_where(mask, quiver_pos_y)
mu = np.ma.masked_where(mask, mu)
mv = np.ma.masked_where(mask, mv)
# # colormapped vectors:
# col = np.sqrt(u**2 + v**2)
# Q = ax.quiver(cell_centers_x[skip], cell_centers_y[skip], u[skip], v[skip], col[skip], units='width', cmap=plt.get_cmap('plasma'), alpha=0.8)
# uniform color:
Q = ax.quiver(mx, my, mu, mv, units='width', color=quivercolor, alpha=0.66)
# qk = plt.quiverkey(Q, 0.9, 0.95, 2, r'$2 \frac{m}{s}$',
# labelpos='E',
# coordinates='figure',
# fontproperties={'weight': 'bold'})
# cax_streams = divider.append_axes("right", size="3%", pad=0.7)
# cbar_streams = plt.colorbar(Q, cax=cax_streams)
# #cbar_streams.ax.set_ylabel('Velocity', rotation=270)
# cbar_streams.set_label('velocity', labelpad=5)
filename_base = basename(filename)
filename_base = os.path.splitext(filename_base)[0]
#filename_base = os.path.splitext(filename)[0]
if not os.path.exists('plots'):
print('--- Make directory...')
os.makedirs('plots', exist_ok=True)
output_dir = 'plots/{}'.format(output_folder)
if not os.path.exists(output_dir):
print('--- Make directory...')
os.makedirs(output_dir, exist_ok=True)
# if not os.path.exists('plots/png'):
# print('--- Make directory...')
# os.makedirs('plots/png', exist_ok=True)
# if not os.path.exists('plots/pdf'):
# print('--- Make directory...')
# os.makedirs('plots/pdf', exist_ok=True)
print('--- Save plot...')
plt.tight_layout()
# get index from filename
# plot_index = map(int, re.findall('\d+', filename_base))
# plot_index = plot_index[0]
# plt.savefig('plots/jpg/frame_%09d.jpg' % plot_index,
# bbox_inches='tight', dpi=200) # dpi=300 200
plt.savefig(os.path.join(output_dir, '{}.jpg'.format(filename_base)),
bbox_inches='tight', dpi=300) # dpi=300 200
# regex = re.compile(r'\d+')
# file_number = regex.search(filename).group(0)
# # regex.findall(filename_base)
# plt.savefig(os.path.join(output_dir, '{}_{}.jpg'.format(scalar_variable, file_number)),
# bbox_inches='tight', dpi=300) # dpi=300 200
#print('Save pdf...')
#plt.savefig('pdfs/'+filename_base+'.pdf', bbox_inches='tight', dpi=300)
# plt.show()
# pp = PdfPages('plots/pdf/' + filename_base + '.pdf')
# pp.savefig()
# pp.close()
plt.close()
print("--- Done!")
return 1
############################################
if __name__ == "__main__":
fn = './data/global_mesh000000.vtu'
# fn = './data/domain_marker000000.vtu'
output_folder = 'marker'
plot_vtk_file(fn, scalar_variable='f', output_folder=output_folder)
# parser = OptionParser()
# parser.add_option("-f", "--file", dest="filename",
# default=False, help="open FILE", metavar="FILE")
# parser.add_option("-q", "--quiet",
# action="store_false", dest="verbose", default=False,
# help="don't print status messages to stdout")
#
# (options, args) = parser.parse_args()
#
# # If no option is given
# if (options.filename == False):
# # If no option is given, plot the whole
# # content of the input-directory:
# input_directory = '../build/bin/out'
# # input_directory = '../build_bc_fix/bin/out'
# # input_directory = '../build_2/bin/out'
# # input_directory = '../build_3/bin/out'
# print("Input directory: %s" % input_directory)
#
# filenames = sorted(glob(input_directory + '/solution*.vt*'))
# interface_filenames = sorted(glob(input_directory + '/interface*.vt*'))
# print("Found the following files:")
# print(filenames)
#
# cmap_binedges = None
# datarange = None
#
# # output = 'velocity'
# # output = 'temperature'
# output = 'density'
# plot_velocity = True
# plot_velocity = False
# # plot_schlieren = True
#
# cmap = plt.cm.Spectral_r
# # cmap = plt.cm.cividis
# # cmap = plt.cm.RdBu_r
# # cmap = plt.cm.Blues
#
# print("get data range...")
# datarange = get_data_range(filenames, output)
# # datarange = get_data_range(filenames, output, separat_ranges=True)
# # cmap_binedges = get_data_bins(filenames, output)
#
# print("Start plotting...")
# output_interval = 50
# output_interval = 1
# # for fn in filenames[::output_interval]:
# for fn, fni in zip(filenames[::output_interval], interface_filenames[::output_interval]):
# print("Plot file {}, {}".format(fn, fni))
# plot_vtk_file(fn, interface_filename=fni, scalar_variable=output, cmap_binedges=cmap_binedges, datarange=datarange, plot_velocity=plot_velocity, cmap=cmap)
# # schlieren_plot_vtk_file(fn, interface_filename=fni, scalar_variable=output, cmap_binedges=cmap_binedges, datarange=datarange, plot_velocity=plot_velocity, cmap=cmap)
#
# video_length = 8.0
# nframes = len(filenames[::output_interval])
# framerate = nframes / video_length
# # framerate = 60
# print('framerate = {}'.format(framerate))
#
# command = ' ffmpeg -framerate {} -pattern_type glob -i "./plots/{}/solution_*.jpg" -filter:v "crop=2*floor(in_w/2):2*floor(in_h/2)" -pix_fmt yuv420p -vcodec h264 -r 60 ./plots/{}.mp4 '.format(framerate, output, output, output)
# # command = ' ffmpeg -framerate {} -pattern_type glob -i "./plots/{}/*.jpg" -filter:v "crop=2*floor(in_w/2):2*floor(in_h/2)" -pix_fmt yuv420p -vcodec h264 -r 60 ./plots/{}.mp4 '.format(framerate, output, output)
# subprocess.call(command, shell=True)
#
# else:
# # If the file is given as a option, plot it
# print("Plot file %s" % options.filename)
# plot_vtk_file(options.filename)
from glob import glob
import os
import sys
# import fileinput
# import re
import shutil
import numpy as np
# import math
# import random
# import matplotlib.pyplot as plt
# from mpl_toolkits.mplot3d import Axes3D
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
# from vtk import *
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from vtk import vtkUnstructuredGridReader
from vtk import vtkXMLUnstructuredGridReader
os.environ["PV_ALLOW_BATCH_INTERACTION"] = "1"
# ------------------------------------------------------------------
def cm2inch(*tupl):
# Convert centimeters to inch
inch = 2.54
if isinstance(tupl[0], tuple):
return tuple(i / inch for i in tupl[0])
else:
return tuple(i / inch for i in tupl)
# ------------------------------------------------------------------
def plot3d_vtu(filename, interface_filename=None, plotname=None, cmap=cm.Spectral_r):
scalar_variable = 'phase'
# interfacecolor = '#0F75FF'
interfacecolor = '#0050CC'
phasecolor = '#73a0e6'
# cmap = sns.light_palette(interfacecolor, as_cmap=True, reverse=False, n_colors=6)
# cmap = cmap(np.linspace(0.0, 0.2, 10))
# cmap = LinearSegmentedColormap.from_list('truncated', cmap)
# cmap = matplotlib.colors.ListedColormap(['white', phasecolor])
# alpha=0.3
# black_cmap = matplotlib.colors.ListedColormap(['white', 'black'])
cmap = cm.get_cmap('Spectral_r')
# Get the filetype
extension = os.path.splitext(filename)[1]
if (extension == '.vtu'):
reader = vtkXMLUnstructuredGridReader()
elif extension == '.vtk':
reader = vtkUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
if interface_filename is not None:
# Get the filetype
interface_extension = os.path.splitext(interface_filename)[1]
if (interface_extension == '.vtu'):
interface_reader = vtkXMLUnstructuredGridReader()
elif interface_extension == '.vtk':
interface_reader = vtkUnstructuredGridReader()
interface_reader.SetFileName(interface_filename)
interface_reader.Update()
interface_output = interface_reader.GetOutput()
print('--- Read cell values...')
n_cell_arrays = reader.GetOutput().GetCellData().GetNumberOfArrays()
n_point_arrays = reader.GetOutput().GetPointData().GetNumberOfArrays()
# time_field = reader.GetOutput().GetFieldData().GetArray(0)
# simulation_time = vtk_to_numpy(time_field)[0]
modified_idx = -1
for i in range(n_cell_arrays):
var_name = reader.GetOutput().GetCellData().GetArrayName(i)
print(var_name)
if (var_name == scalar_variable):
scalar_idx = i
if (var_name == 'modified'):
modified_idx = i
# scalar_vtk_array = reader.GetOutput().GetCellData().GetArray(scalar_idx)
# reader.GetOutput().GetCellData().SetActiveScalars(scalar_variable)
# reader.GetOutput().GetCellData().SetActiveScalars('phase')
# output = interface_reader.GetOutput()
output = reader.GetOutput()
bounds = reader.GetOutput().GetBounds()
center = reader.GetOutput().GetCenter()
scalar_range = output.GetScalarRange()
print(scalar_range)
# ----------------------------------------
# cmap_inst.set_bad(color='white')
vmin, vmax = scalar_range
lut = vtk.vtkLookupTable()
lut.SetTableRange(vmin, vmax)
# lut.SetNanColor(.1, .5, .99, 1.0)
n = 256
lut.SetNumberOfTableValues(n)
lut.Build()
# for i in range(n):
for i, val in enumerate(np.linspace(0, 1, n)):
# R,G,B = colors.colorMap(i, 'jet', 0, n)
rgb_val = cmap(val)
lut.SetTableValue(i, rgb_val[0], rgb_val[1], rgb_val[2], 1)
# ----------------------------------------
mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(output)
# mapper.SetInputConnection(output_port)
mapper.SetScalarRange(scalar_range)
mapper.SetScalarModeToUseCellData()
# mapper.SelectColorArray('phase')
mapper.SetLookupTable(lut)
# ----------------------------------------
#
# interface_mapper = vtk.vtkDataSetMapper()
# interface_mapper.SetInputData(interface_output)
# interface_actor = vtk.vtkActor()
# interface_actor.SetMapper(interface_mapper)
# interface_actor.GetProperty().EdgeVisibilityOn()
# interface_actor.GetProperty().SetLineWidth(1.0)
# interface_actor.GetProperty().SetColor(matplotlib.colors.to_rgb(interfacecolor))
# # interface_actor.GetProperty().SetSpecular(0)
# # interface_actor.GetProperty().SetSpecularPower(0)
# interface_actor.GetProperty().SetOpacity(0.2)
# interface_actor.GetProperty().SetRepresentationToWireframe()
# # interface_actor.GetProperty().SetAmbient(0.1)
# # interface_actor.GetProperty().SetDiffuse(1)
# ----------------------------------------
#
# plane = vtk.vtkPlane()
# plane.SetOrigin(output.GetCenter())
# plane.SetNormal(0, 0, 1)
# plane_cut = vtk.vtkCutter()
# plane_cut.GenerateTrianglesOff()
# plane_cut.SetInputData(output)
# plane_cut.SetCutFunction(plane)
# cut_mapper = vtk.vtkDataSetMapper()
# cut_mapper.SetInputConnection(plane_cut.GetOutputPort())
# cut_mapper.SetScalarRange(scalar_range)
# cut_mapper.ScalarVisibilityOn()
# cut_mapper.SetScalarModeToUseCellData()
# cut_mapper.SetLookupTable(lut)
# ----------------------------------------
#
# extracter = vtk.vtkExtractGeometry()
#
# extracter.SetInputData(output);
# extracter.SetImplicitFunction(plane);
# extracter.ExtractInsideOn();
# extracter.ExtractBoundaryCellsOn();
#
# extract_mapper = vtk.vtkDataSetMapper()
# extract_mapper.SetInputConnection(extracter.GetOutputPort())
# extract_mapper.SetScalarRange(scalar_range)
# extract_mapper.ScalarVisibilityOn()
# extract_mapper.SetScalarModeToUseCellData()
# extract_mapper.SetLookupTable(lut)
# # ----------------------------------------
#
# clipper = vtk.vtkClipDataSet()
# clipper.SetClipFunction(plane)
# clipper.SetInputData(reader.GetOutput())
# clipper.InsideOutOn()
# clipper.SetValue(0.0)
# clipper.GenerateClippedOutputOn()
# clipper.Update()
#
# clip_mapper = vtk.vtkDataSetMapper()
# clip_mapper.SetInputData(clipper.GetOutput())
#
# # ----------------------------------------
main_mapper = mapper
# Create the Actor
actor = vtk.vtkActor()
actor.SetMapper(main_mapper)
# actor.SetMapper(clip_mapper)
# actor.SetMapper(extract_mapper)
# actor.SetMapper(cut_mapper)
actor.GetProperty().EdgeVisibilityOn()
actor.GetProperty().SetLineWidth(1.0)
actor.GetProperty().SetSpecular(0)
actor.GetProperty().SetSpecularPower(0)
# actor.GetProperty().SetOpacity(0.5)
# actor.GetProperty().SetRepresentationToWireframe()
actor.GetProperty().SetAmbient(1)
actor.GetProperty().SetDiffuse(0)
# ----------------------------------------
colorbar = vtk.vtkScalarBarActor()
colorbar.SetLookupTable(main_mapper.GetLookupTable())
colorbar.SetTitle(scalar_variable)
# colorbar.SetOrientationToHorizontal()
colorbar.SetOrientationToVertical()
# colorbar.GetLabelTextProperty().SetColor(0,0,0)
# colorbar.GetTitleTextProperty().SetColor(0,0,0)
tprop = vtk.vtkTextProperty()
# tprop.SetFontFamily(vtk.VTK_FONT_FILE)
# # tprop.SetFontFamilyAsString('Linux Biolinum O')
# # tprop.SetFontFile('./LibertinusSans-Bold.otf')
# tprop.SetFontFile('./LibertinusSans-Regular.otf')
tprop.SetColor(0,0,0)
# tprop.BoldOn()
colorbar.SetLabelTextProperty(tprop)
colorbar.SetTitleTextProperty(tprop)
# position it in window
coord = colorbar.GetPositionCoordinate()
coord.SetCoordinateSystemToNormalizedViewport()
coord.SetValue(0.9,0.2)
colorbar.SetWidth(.075)
colorbar.SetHeight(.5)
# ----------------------------------------
# backface = vtk.vtkProperty()
colors = vtk.vtkNamedColors()
# backface.SetColor(colors.GetColor3d("tomato"))
# actor.SetBackfaceProperty(backface)
# ----------------------------------------
# camera = renderer.GetActiveCamera()
camera = vtk.vtkCamera()
# camera.SetPosition(0, 0, 10)
# camera.SetPosition(-1, -1, 5)
# camera.SetFocalPoint(0, 0, 0)
camera.SetViewUp(0, 1, 0)
camera.ParallelProjectionOn()
# camera.SetParallelScale(7)
# renderer.ResetCamera()
spacing = [1.1,1.1]
xc = center[0] + 0.5*(bounds[0] + bounds[1])*spacing[0]
yc = center[1] + 0.5*(bounds[2] + bounds[3])*spacing[1]
xd = (bounds[1] - bounds[0] + 1)*spacing[0]
yd = (bounds[3] - bounds[2] + 1)*spacing[1]
d = camera.GetDistance()
camera.SetParallelScale(0.5*yd)
camera.SetFocalPoint(xc,yc,0.0)
camera.SetPosition(xc,yc,+d)
# ----------------------------------------
# renderer = vtk.vtkRenderer()
#
light = vtk.vtkLight()
light.SetFocalPoint(0,0,0)
light.SetPosition(-1, -1, 10)
# ----------------------------------------
axis_actor = vtk.vtkCubeAxesActor2D()
axis_actor.SetBounds(output.GetBounds())
axis_actor.SetLabelFormat("%6.4g")
axis_actor.SetFlyModeToOuterEdges()
axis_actor.SetFontFactor(1.5)
# tprop = vtk.vtkTextProperty()
# tprop.SetFontFamilyAsString('Linux Biolinum O')
# tprop.SetColor(0,0,0)
# tprop.ShadowOn()
axis_actor.ScalingOn()
axis_actor.SetAxisTitleTextProperty(tprop)
axis_actor.SetAxisLabelTextProperty(tprop)
axis_actor.SetFlyModeToOuterEdges()
# axis_actor.SetFlyModeToClosestTriad()
# axis_actor.YAxisVisibilityOn()
axis_actor.YAxisVisibilityOff()
# axis_actor = vtk.vtkCubeAxesActor()
# axis_actor.SetUseTextActor3D(1)
# axis_actor.SetBounds(output.GetBounds())
# axis_actor.SetCamera(renderer.GetActiveCamera())
# axis_actor.GetTitleTextProperty(0).SetColor(0,0,0)
# axis_actor.GetTitleTextProperty(0).SetFontSize(48)
# axis_actor.GetLabelTextProperty(0).SetColor(0,0,0)
# axis_actor.GetTitleTextProperty(1).SetColor(0,0,0)
# axis_actor.GetLabelTextProperty(1).SetColor(0,0,0)
# axis_actor.GetTitleTextProperty(2).SetColor(0,0,0)
# axis_actor.GetLabelTextProperty(2).SetColor(0,0,0)
# axis_actor.DrawXGridlinesOn()
# axis_actor.DrawYGridlinesOn()
# axis_actor.DrawZGridlinesOn()
# axis_actor.SetGridLineLocation(axis_actor.VTK_GRID_LINES_FURTHEST)
# axis_actor.XAxisMinorTickVisibilityOff()
# axis_actor.YAxisMinorTickVisibilityOff()
# axis_actor.ZAxisMinorTickVisibilityOff()
# axis_actor.SetFlyModeToStaticEdges()
# ----------------------------------------
# # Create a text actor.
# txt_actor = vtk.vtkTextActor()
# txt_actor.SetInput('t = {:.1e}'.format(simulation_time))
# txt_actor.SetTextProperty(tprop)
# txt_actor.GetTextProperty().SetFontSize(16)
# # txt_actor.SetTextScaleModeToProp()
# txt_actor.SetTextScaleModeToViewport()
# # txt_actor.SetTextScaleModeToNone()
# txt_actor.GetTextProperty().SetJustificationToRight()
# txt_actor.GetTextProperty().SetVerticalJustificationToTop()
# txt_actor.GetTextProperty().UseTightBoundingBoxOn()
#
# # txt_actor.SetDisplayPosition(20, 30)
# coord = txt_actor.GetPositionCoordinate()
# coord.SetCoordinateSystemToNormalizedViewport()
# coord.SetValue(0.975,0.975)
# # txt_actor.SetWidth(1)
# # txt_actor.SetHeight(1)
# ----------------------------------------
renderer = vtk.vtkRenderer()
renderer.SetActiveCamera(camera)
# renderer.GetActiveCamera().Zoom(3.33)
# renderer.GetActiveCamera().Dolly(1.0)
renderer.AddActor(actor)
# renderer.AddActor(colorbar)
axis_actor.SetCamera(renderer.GetActiveCamera())
renderer.AddActor(axis_actor)
# renderer.AddActor(txt_actor)
renderer.AddLight(light)
# renderer.AddActor(interface_actor)
renderer.SetBackground(1, 1, 1)
# renderer.SetBackground(colors.GetColor3d("Wheat"))
renderer.UseHiddenLineRemovalOn()
# renderer.ResetCameraClippingRange()
# renderer.SetUseDepthPeeling(1)
# renderer.SetOcclusionRatio(0.1)
# renderer.SetMaximumNumberOfPeels(100)
# render_window.SetMultiSamples(0)
# render_window.SetAlphaBitPlanes(1)
# Create the RendererWindow
renderer_window = vtk.vtkRenderWindow()
renderer_window.SetOffScreenRendering(1)
renderer_window.AddRenderer(renderer)
# renderer_window.Render()
# create a renderwindowinteractor
renderer_window.Render()
# renderer_window.SetWindowName('Plot')
# renderer_window.SetSize(600, 300)
# renderer_window.SetSize(500, 500)
basesize = 1000
renderer_window.SetSize(int(1.1*basesize), basesize)
# renderer_window.SetSize(2000, 2000)
# win_size = renderer_window.GetSize()
# factor = 2
# renderer_window.SetSize(factor * win_size[0], factor * win_size[1])
# ----------------------------------------
# write plot
renderer_window.Render()
windowto_image_filter = vtk.vtkWindowToImageFilter()
windowto_image_filter.SetInput(renderer_window)
rgba = False
if rgba:
windowto_image_filter.SetInputBufferTypeToRGBA()
windowto_image_filter.Update()
else:
windowto_image_filter.SetInputBufferTypeToRGB()
windowto_image_filter.ReadFrontBufferOff()
windowto_image_filter.Update()
# ----------------------------------------
filename = plotname + '.jpg'
writer = vtk.vtkJPEGWriter()
# writer = vtk.vtkPNGWriter()
writer.SetFileName(filename)
writer.SetInputConnection(windowto_image_filter.GetOutputPort())
writer.Write()
return 1
# Create the RendererWindowInteractor and display the vtk_file
# interactor = vtk.vtkRenderWindowInteractor()
# interactor.SetRenderWindow(renderer_window)
interactor.Initialize()
interactor.Start()
print("--- Done!")
return 1
# ------------------------------------------------------------------------------
if __name__ == "__main__":
fn = './data/global_mesh000000.vtu'
fn = './data/domain_marker000000.vtu'
plot3d_vtu(fn, interface_filename=None, plotname='./plots/vtu_plot')
# basepath = '../build/bin'
# # basepath = '../build_benchmarks/bin'
# # basepath = '.'
# paths = os.listdir(basepath)
# # print(paths)
# folders = [os.path.basename(p) for p in paths if os.path.isdir(
# os.path.join(basepath, p))]
# # print(folders)
# for idx, f in enumerate(folders):
# print("{}: {}".format(idx, f))
# idx = int(input("Choose dataset directory: "))
# # idx = 1
# main_dir = folders[idx]
# benchmark_name = folders[idx]
# # main_dir = os.path.join('.', folders[idx])
# print("Directory: {}".format(main_dir))
# main_dir = os.path.join(basepath, main_dir)
# print("Directory: {}".format(main_dir))
#
# all_paths = []
# all_paths.extend(glob(main_dir + '/benchmark*'))
# paths = [p for p in all_paths if os.path.isdir(p)]
# result = []
# for p in paths:
# if any(fname.startswith('statistics') for fname in os.listdir(p)):
# result.append(p)
# paths = result
# paths = sorted(paths)
#
# # print(paths)
# npaths = len(paths)
# print('There are %d benchmark folders.' % npaths)
# # sys.exit()
#
#
# # ==========================================================================
# # plot grids
#
# paths = [p for p in all_paths if os.path.isdir(p)]
# paths = sorted(paths)
# for idx, f in enumerate(paths):
# print("{}: {}".format(idx, os.path.basename(f) ))
# idx = int(input("Choose data to plot: "))
# # idx = 0
# plot_path = paths[idx]
# input_directory = os.path.join(plot_path, 'vtk')
#
# print("Input directory: %s" % input_directory)
#
# filenames = sorted(glob(input_directory + '/solution*.vt*'))
# interface_filenames = sorted(glob(input_directory + '/interface*.vt*'))
# print("Found the following files:")
# print(filenames)
#
# cmap = plt.cm.Spectral_r
# # cmap = plt.cm.cividis
# # cmap = plt.cm.RdBu_r
# # cmap = plt.cm.Blues
#
# plot_output_dir = os.path.join(main_dir, 'plots')
# os.makedirs(plot_output_dir, exist_ok=True)
# shutil.copy(plot_path + '/log.json', plot_output_dir + '/log.json')
#
# print("Start plotting...")
# output_interval = 1
# # output_interval = 50
#
# file_list = list(zip(filenames[::output_interval], interface_filenames[::output_interval]))
#
# np.random.shuffle(file_list)
#
# # for fn in filenames[::output_interval]:
# # for fn, fni in zip(filenames[::output_interval], interface_filenames[::output_interval]):
#
# # file_list = [random.sample(file_list, k=1)]
#
# for fn, fni in file_list:
# print("Plot file {}, {}".format(fn, fni))
# plotname = os.path.basename(fn)
# plotname = os.path.splitext(plotname)[0]
# plotname = os.path.join(plot_output_dir, plotname)
# # if os.path.isfile(plotname):
# # continue
# print(plotname)
# plot3d_vtu(fn, interface_filename=fni, plotname=plotname, cmap=cmap)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment