setup new TPR layered soil examples

Plotskripte_Jim/links

+https://lorensen.github.io/VTKExamples/site/Python/

Plotskripte_Jim/plot2d_vtu.py

+#!/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)

Plotskripte_Jim/plot3d_vtu.py

+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)

Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-different-permeabilties.py

@@ -54,10 +54,10 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = {0: 0.0}
-# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
+# starttimes = {0: 0.0}
+starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
 timestep_size = 0.001
-number_of_timesteps = 1500
+number_of_timesteps = 1
 
 # LDD scheme parameters  ######################################################
 Lw1 = 0.01  # /timestep_size
@@ -69,41 +69,41 @@ Lnw2 = Lw2
 Lw3 = 0.002 # /timestep_size
 Lnw3 = 0.002
 
-Lw4 = 0.002 # /timestep_size
-Lnw4 = 0.002
+Lw4 = 0.0001 # /timestep_size
+Lnw4 = 0.0001
 
-Lw5 = 0.002 # /timestep_size
-Lnw5 = 0.002
+Lw5 = 0.0001 # /timestep_size
+Lnw5 = 0.0001
 
-Lw6 = 0.002  # /timestep_size
-Lnw6 = 0.002
+Lw6 = 0.0001  # /timestep_size
+Lnw6 = 0.0001
 
 lambda12_w = 1
-lambda12_nw = 0.4
+lambda12_nw = 1
 
 lambda23_w = 1
-lambda23_nw = 0.4
+lambda23_nw = 1
 
 lambda24_w = 1
-lambda24_nw= 0.4
+lambda24_nw= 1
 
 lambda25_w= 1
-lambda25_nw= 0.4
+lambda25_nw= 1
 
 lambda34_w = 1
-lambda34_nw = 0.4
+lambda34_nw = 1
 
 lambda36_w = 1
-lambda36_nw = 0.4
+lambda36_nw = 1
 
 lambda45_w = 1
-lambda45_nw = 0.4
+lambda45_nw = 1
 
 lambda46_w = 1
-lambda46_nw = 0.4
+lambda46_nw = 1
 
 lambda56_w = 1
-lambda56_nw = 0.4
+lambda56_nw = 1
 
 include_gravity = False
 debugflag = False
@@ -113,10 +113,10 @@ analyse_condition = False
 # when number_of_timesteps is high, it might take a long time to write all
 # timesteps to disk. Therefore, you can choose to only write data of every
 # plot_timestep_every timestep to disk.
-plot_timestep_every = 3
+plot_timestep_every = 1
 # Decide how many timesteps you want analysed. Analysed means, that
 # subsequent errors of the L-iteration within the timestep are written out.
-number_of_timesteps_to_analyse = 5
+number_of_timesteps_to_analyse = 1
 
 # fine grained control over data to be written to disk in the mesh study case
 # as well as for a regular simuation for a fixed grid.
@@ -225,9 +225,9 @@ porosity = {
     1: 0.57,  #0.2,  # Clayey gravels, clayey sandy gravels
     2: 0.57,  #0.22, # Silty gravels, silty sandy gravels
     3: 0.005,  #0.57, # Clayey sands
-    4: 0.005,  #0.2 # Silty or sandy clay
-    5: 0.005,  #
-    6: 0.005,  #
+    4: 0.0005,  #0.2 # Silty or sandy clay
+    5: 0.0001,  #
+    6: 0.001,  #
 }
 
 # subdom_num : subdomain L for L-scheme
@@ -292,8 +292,8 @@ intrinsic_permeability = {
     1: 0.1,  # sand
     2: 0.1,  # sand, there is a range
     3: 0.001,  #10e-2,  # clay has a range
-    4: 0.001,  #10e-3
-    5: 0.001,  #10e-2,  # clay has a range
+    4: 0.0005,  #10e-3
+    5: 0.0001,  #10e-2,  # clay has a range
     6: 0.001,  #10e-3
 }
 

Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py

@@ -61,7 +61,6 @@ timestep_size = 0.001
 number_of_timesteps = 1000
 
 # LDD scheme parameters  ######################################################
-
 Lw1 = 0.01  # /timestep_size
 Lnw1 = Lw1