output binary and add layer mask

config.yaml

 filename_in: "/home/jarunanp/work/geg/19_04_Justin_Paraview/Production_openBC/Results/openBC_prod_10yr.csv"
-filename_out: "/home/jarunanp/work/geg/19_04_Justin_Paraview/output/openBC_prod_10yr.vtk"
+filename_out: "/home/jarunanp/work/geg/19_04_Justin_Paraview/output/openBC_prod_10yr"
+
+# dimension
 npoint_d1: 80
 npoint_d2: 40
 npoint_d3: 15
+
+# layer thickness only in z direction
+layer_top: 5
+layer_middle: 5
+layer_bottom: 5
+

results_to_vtk.py

@@ -2,20 +2,35 @@ import os
 import sys
 import numpy as np
 import yaml
-
+from pyevtk.hl import gridToVTK 
+import random as rnd
 
 def read_config_file(filename):
     """ read input values """
 
-    print(" ... reading an input file", filename)
+    print("... reading an input file", filename)
     with open(filename, "r") as f:
         config = yaml.load(f)
 
     return config
 
 
-def read_results(config):
-    """ read mesh and results """
+def create_layer_mask(npoint_d1, npoint_d2, npoint_d3, top, middle, bottom):
+    layer_mask = np.zeros((npoint_d1, npoint_d2, npoint_d3))
+    for k in range(npoint_d3): 
+        for j in range(npoint_d2):
+            for i in range(npoint_d1):
+                if (k < bottom):
+                    layer_mask[i,j,k] = 0
+                elif (k < middle+bottom):
+                    layer_mask[i,j,k] = 2
+                elif (k < top+middle+bottom):
+                    layer_mask[i,j,k] = 4
+    return layer_mask
+
+
+def output_ascii_vtk(config):
+    """ read mesh and results and output in ASCII format """
 
     print(" ... reading a result file", config["filename_in"])
     # Open the file and read data
@@ -45,8 +60,10 @@ def read_results(config):
     ncell_d2 = npoint_d2-1
     ncell_d3 = npoint_d3-1
     ncells = ncell_d1*ncell_d2*ncell_d3
+
+    filename_out = "{}.vtk".format(config["filename_out"])
     # Write output
-    with open(config["filename_out"], "w") as g:
+    with open(filename_out, "w") as g:
         g.write("# vtk DataFile Version 2.0\n")
         g.write("Mesh data\n")
         g.write("ASCII\n")
@@ -76,10 +93,73 @@ def read_results(config):
                g.write("{} ".format(value))
             g.write("\n")
 
+
+def output_binary_vtk(config):
+    """ read mesh and results and output in BINARY format """
+    print("... reading a result file", config["filename_in"])
+
+    npoint_d1 = config["npoint_d1"]
+    npoint_d2 = config["npoint_d2"]
+    npoint_d3 = config["npoint_d3"]
+    npoints = npoint_d1*npoint_d2*npoint_d3
+    #
+    ncell_d1 = npoint_d1-1
+    ncell_d2 = npoint_d2-1
+    ncell_d3 = npoint_d3-1
+    ncells = ncell_d1*ncell_d2*ncell_d3
+
+    # Open the file and read data
+    with open(config["filename_in"], "r") as f:
+        data=f.readlines()
+
+    header = data.pop(0)
+    header = header.replace('\"','')
+    header = header.replace(' ','')
+    header = header.strip().split(",")
+    header.pop(0)
+    header.pop(0)
+    header.pop(0)
+    #
+    nresults = len(header)
+    print("... reading {:d} results: ".format(nresults))
+    for i in range(nresults):
+       print("   {}".format(header[i]))
+    #
+    results = {}
+    for i in range(nresults):
+        key = header[i]
+        results[key] = np.zeros((npoint_d1, npoint_d2, npoint_d3))
+
+    # initialize coordinates
+    x = np.zeros((npoint_d1, npoint_d2, npoint_d3)) 
+    y = np.zeros((npoint_d1, npoint_d2, npoint_d3)) 
+    z = np.zeros((npoint_d1, npoint_d2, npoint_d3)) 
+
+    # Write output
+    for k in range(npoint_d3): 
+        for j in range(npoint_d2):
+            for i in range(npoint_d1):
+               ii = i + j*npoint_d1 + k*npoint_d1*npoint_d2
+               line = data[ii]
+               line = line.replace('NAN','0.0')
+               values = line.strip().split(",")
+               x[i,j,k] = float(values[0])
+               y[i,j,k] = float(values[1])
+               z[i,j,k] = float(values[2])
+               # store results to arrays
+               for l in range(nresults):
+                   key = header[l]
+                   results[key][i,j,k] = float(values[l+3])
+
+    layer_mask = create_layer_mask(npoint_d1, npoint_d2, npoint_d3, config["layer_top"], config["layer_middle"], config["layer_bottom"])
+    results['layer_mask'] = layer_mask
+
+    gridToVTK(config["filename_out"], x, y, z, pointData = results)
+
 if __name__ == "__main__":
 
-    # read config file
+   # read config file
     config = read_config_file("config.yaml")
     
     # read result file
-    read_results(config)
+    output_binary_vtk(config)