From ecc5d3f5665340370217e71ba0ab4a85d86e1df7 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Mon, 12 Oct 2020 10:47:47 +0200
Subject: [PATCH] remove Jim's python plotscript due to copyright concerns
---
Meshes/Plotscript/links | 1 -
Meshes/Plotscript/plot2d_vtu.py | 502 -------------------------------
Meshes/Plotscript/plot3d_vtu.py | 510 --------------------------------
3 files changed, 1013 deletions(-)
delete mode 100644 Meshes/Plotscript/links
delete mode 100644 Meshes/Plotscript/plot2d_vtu.py
delete mode 100644 Meshes/Plotscript/plot3d_vtu.py
diff --git a/Meshes/Plotscript/links b/Meshes/Plotscript/links
deleted file mode 100644
index bffbaf4..0000000
--- a/Meshes/Plotscript/links
+++ /dev/null
@@ -1 +0,0 @@
-https://lorensen.github.io/VTKExamples/site/Python/
diff --git a/Meshes/Plotscript/plot2d_vtu.py b/Meshes/Plotscript/plot2d_vtu.py
deleted file mode 100644
index 08d5b09..0000000
--- a/Meshes/Plotscript/plot2d_vtu.py
+++ /dev/null
@@ -1,502 +0,0 @@
-#!/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)
diff --git a/Meshes/Plotscript/plot3d_vtu.py b/Meshes/Plotscript/plot3d_vtu.py
deleted file mode 100644
index 683a186..0000000
--- a/Meshes/Plotscript/plot3d_vtu.py
+++ /dev/null
@@ -1,510 +0,0 @@
-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)
--
GitLab