Improvements to plotting function for flow-over-heated-plate

flow-over-heated-plate/plot-final-interface-temperature-requirements.txt

@@ -1,3 +1,4 @@
 vtk
 numpy
 matplotlib
+os

flow-over-heated-plate/plot-final-interface-temperature.py

@@ -1,68 +1,65 @@
+#!/usr/bin/env python3
 import vtk
 from matplotlib import pyplot as plt
 import numpy as np
+import os
 
 
 def vtk_to_dict(case):
-    vtkFileName = "solid-{}/precice-exports/Fluid-Mesh-Solid.dt100.vtk".format(case)
+    vtkFileName = "solid-{}/precice-exports/Fluid-Mesh-Solid.dt100.vtu".format(
+        case)
+    if not os.path.exists(vtkFileName):
+        print("No file found for " + vtkFileName)
+        return {}  # return empty dict if file not found
+
     # read the vtk file as an unstructured grid
-    reader = vtk.vtkUnstructuredGridReader()
+    reader = vtk.vtkXMLUnstructuredGridReader()
     reader.SetFileName(vtkFileName)
-    reader.ReadAllVectorsOn()
-    reader.ReadAllScalarsOn()
     reader.Update()
 
     # obtain the data
     data = reader.GetOutput()
     n_data = data.GetPointData().GetNumberOfTuples()
 
-    name = "Temperature"
-    data_names = []
-    i = 0
-    max_i = data.GetPointData().GetNumberOfArrays()
-    while i < max_i:
-        this_data_name = data.GetPointData().GetArray(i).GetName()
-        data_names.append(this_data_name)
-        if (this_data_name == name):
-            data_id = i
-            break
-        i += 1
-
     data_dict = {}
 
-    if not data_id:
-        raise Exception(
-            "For file {} name {} not found. Only the following names are available: {}. "
-            "Aborting!".format(vtkFileName, name, data_names))
     for i in range(n_data):
-        data_dict[data.GetPoint(i)] = data.GetPointData().GetArray(data_id).GetValue(i)
-
+        data_dict[data.GetPoint(i)] = data.GetPointData().GetArray(
+            "Temperature").GetValue(i)
     return data_dict
 
 
-cases = []
-cases.append('fenics')
-cases.append('openfoam')
-cases.append('nutils')
-cases.append('dunefem')
+def main():
+    case_labels = {
+        'fenics': 'Fluid-FEniCS',
+        'openfoam': 'Fluid-OpenFOAM',
+        'nutils': 'Fluid-Nutils',
+        'dunefem': 'Fluid-DuneFem'}
+    styles = [':', '-', '--']
+    colors = ['r', 'b', 'g', 'k']
+
+    for i, case in enumerate(case_labels.keys()):
+        case_data = vtk_to_dict(case)
+        if not case_data:
+            continue
+        x, t = [p[0] for p in case_data.keys()], np.array(
+            list(case_data.values()))
+
+        # sort by x
+        combined = sorted(zip(x, t))
+        x, t = zip(*combined)
+        x = np.array(x)
+        t = np.array(t)
+
+        theta = (t - 300) / (310 - 300)
+        plt.plot(x, theta, colors[i % 4] + styles[i %
+                 3], label=case_labels[case])
 
-case_labels = {
-    'fenics': 'OpenFOAM-FEniCS',
-    'openfoam': 'OpenFOAM-OpenFOAM',
-    'nutils': 'OpenFOAM-Nutils',
-    'dunefem': 'OpenFOAM-DuneFem'}
-styles = [':', '-', '--']
-colors = ['r', 'b', 'g', 'k']
-i = 0
+    plt.ylabel("Theta")
+    plt.xlabel("x-coordinate along coupling interface")
+    plt.legend()
+    plt.show()
 
-for case in cases:
-    case_data = vtk_to_dict(case)
-    x, t = [p[0] for p in case_data.keys()], np.array(list(case_data.values()))
-    theta = (t - 300) / (310 - 300)
-    plt.plot(x, theta, colors[i % 4] + styles[i % 3], label=case_labels[case])
-    i += 1
 
-plt.ylabel("Theta")
-plt.xlabel("x-coordinate along coupling interface")
-plt.legend()
-plt.show()
+if __name__ == '__main__':
+    main()