diff --git a/Meshes/Plotscript/links b/Meshes/Plotscript/links deleted file mode 100644 index bffbaf4e36b0d4ad55f34b914c8760e1fbc1909b..0000000000000000000000000000000000000000 --- 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 08d5b09c1ed72ff38e36a88092545e930c73eaf5..0000000000000000000000000000000000000000 --- 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 683a186c9854b7d57638e1ce481e9302600c8189..0000000000000000000000000000000000000000 --- 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)