CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   Postprocessing .vtk file OpenFOAM (https://www.cfd-online.com/Forums/openfoam/223667-postprocessing-vtk-file-openfoam.html)

kaboka January 21, 2020 12:06

Postprocessing .vtk file OpenFOAM
 
Hello everyone,

I am simulating a sloshing case inside the 3D cylindirical tank. I used the following code in my control dict file and got the free surface points as a .vtk file. However I have these .vtk files for my every timesteps for all the points of free surface. What I want is to plot a graph of the deflection of the point which is at the leftmost corner for all time steps. Could anybody help me how to post process this point for every timestep and plot that graph in python or excel please ? or is there any other method to do it?

Thank in advance
HTML Code:

functions
{
  freeSurface
  { 
      type            surfaces;
      functionObjectLibs
      ( 
          "libsampling.so"
      ); 
      outputControl  outputTime;
      outputInterval  1; 
      surfaceFormat  vtk;
      fields
      ( 
          alpha.water
      ); 
      surfaces
      ( 
          freeSurface
          { 
              type        isoSurfaceCell;
              isoField    alpha.water;
              isoValue    0.5;
              interpolate false;
              regularise  false;
          } 
      ); 
      interpolationScheme cell;
  } 
);


Antimony January 22, 2020 01:29

Hi,

You should be able to do this using the vtk library/package in python. I had to do something similar and my approach was to use python vtk library to read in and extract information from .vtk files

Hope this helps.

Cheers,
Antimony

kaboka January 22, 2020 10:11

Hi Antimony,

Thanks for the answer. I tried it with python but somehow it doesn't work. If it's possible could you share your python code with me? Thanks a lot in advance. I also put the python code below.
HTML Code:

#!/usr/bin/python

# elevationVsTime
# Read VTK files with isosurface
# Track one (x,z) coordinate in time

import os
import re
from vtk import *
from optparse import OptionParser
from numpy import *

print ("elevationVsTime v0.2")

# Read command line arguments
parser = OptionParser()

parser.add_option("-f","--input-file",dest="coords",type="string",help="Filename containing coordinates",metavar="FILE",default="coords")

(options, args) = parser.parse_args()

print ("Reading coordinates from file \"",options.coords,"\"")

f = loadtxt(options.coords)


#- Search starting point
x = f[:,0]
y = x*0+0.7
z = f[:,1]

points = len(x)

print (points," points found")

#- List input points
for i in range(0, points-1):
 print ("Point ",i," (",x[i],",",z[i],")")

# Import timedirectories
# read the vtk directory and get all the time steps and return list
basedir = "freeSurface/"

timesteps=[]
for root,dir,file in os.walk(basedir,True):
 p,time = os.path.split(root)
 if (bool(re.search("^[0-9.]",time))):
  timesteps.append(time)


filename = file[0]
basename = 'elevationVsTime'
timesteps = sorted(timesteps) # This sorts alphabetically

# "Progress-bar" necessity -
Ntimes = len(timesteps)
frac = round(Ntimes/10.0);
counter = 0

#- Time-info - known bug here
print ("Timerange: [",timesteps[0],",",timesteps[Ntimes-1],"]")

## Read
for ts in timesteps:
  #- Counter info -- this info is approximate
  counter = counter+1
  if ( counter%frac == 0 ):
  print (round(counter/frac)*10,"%")

  #- Read VTK file
  readfile = basedir + ts + "/" + filename
  reader = vtkPolyDataReader()
  reader.SetFileName(readfile)
  reader.Update()
  output = reader.GetOutput()

  #- For each timestep, find the closest coordinate on the 0/ plane
  for i in range(points):
  writefile = basename + str(i)
  file = open(writefile,'a+')

  #- Coordinate to find closest point
  xfind = [x[i], y[i], z[i]]

  #- Find point
  p = output.FindPoint(xfind)
 
  #- Coordinate of point
  xfound = output.GetPoint(p)
  y[i] = xfound[1]

  #- Write to file
  print >> file, ts, y[i]
  file.close()


kaboka March 24, 2020 04:56

Quote:

Originally Posted by Antimony (Post 755268)
Hi,

You should be able to do this using the vtk library/package in python. I had to do something similar and my approach was to use python vtk library to read in and extract information from .vtk files

Hope this helps.

Cheers,
Antimony

Could you help me with the programming, please?

Antimony March 31, 2020 23:18

Hi,

What is the issue that you are facing with the code? Can you explain in detail?

Cheers,
Antimony

kaboka April 1, 2020 07:05

Quote:

Originally Posted by Antimony (Post 763744)
Hi,

What is the issue that you are facing with the code? Can you explain in detail?

Cheers,
Antimony

Hi Antimony,

I am getting this error. I think I have to define a coordinate system that is inside the coords.txt but I do not know what the code exactly wants from me.
PHP Code:

OSErrorcoords not found

I got the code from the website and it says that;
Quote:

Gravity is assumed to act in the y-direction
A file 'coords' with two columns (x, z) is expected as probing coordinates
This is the website:
HTML Code:

https://openfoamwiki.net/index.php/Tip_Surface_elevation_in_time
To be honest, I am not good at python or programming, thus I am super confused about that coords thing.

I would be appreciated if you help me. Thanks in advance

zhangyan April 1, 2020 09:11

Hello,
I am really interested in using python to plot with the vtk file input.
Could you please offer an example?


Quote:

Originally Posted by Antimony (Post 755268)
Hi,

You should be able to do this using the vtk library/package in python. I had to do something similar and my approach was to use python vtk library to read in and extract information from .vtk files

Hope this helps.

Cheers,
Antimony


crubio.abujas May 21, 2020 09:45

Possible Code solution?
 
Hi kaboka,

Quote:

Originally Posted by kaboka (Post 763806)
I am getting this error. I think I have to define a coordinate system that is inside the coords.txt but I do not know what the code exactly wants from me.
PHP Code:

OSErrorcoords not found


It seems that the code is looking for a coordinate file and cannot find it. The code is intended to be used in command line fashion, so you shall be executing something similar to:
Code:

elevationVsTime --input-file coordinate_file
where the coordinate_file shall contain a list of points you're interested on.
It assumes that the points are in (X,Z) format, probably because that was the special need of the author. The code use these points to check against the different surfaces what are the closest ones. Then the y coordinate is saved into a file.

Quote:

Originally Posted by kaboka (Post 763806)
I got the code from the website and it says that;

This is the website:
HTML Code:

https://openfoamwiki.net/index.php/Tip_Surface_elevation_in_time
To be honest, I am not good at python or programming, thus I am super confused about that coords thing.

You said that you're interested in the deflection of the leftmost point, maybe
is it possible to simplify the problem to getting the high of the highest
point? I've wrote a simple code for doing so, maybe you can tweak it for your
specific purposes. You may have to change the data on the beginning for your
case. In the core you have a np.array of the points in the surface, so you can
do whatever manipulation you may consider necessary.

Code:

#!/usr/bin/python3

import os, sys, re

import vtk
import matplotlib.pyplot as plt
import numpy as np

### - DATA - ###
BASEDIR = '../postProcessing/freeSurface'
output_file = 'results.out'

print('- - - - - Searching VTK surfaces - - - - -')

# Ensure that the folder exists
if not os.path.exists(BASEDIR):
    print('ERROR: No folder "{}" found'.format(BASEDIR))
    sys.exit(1)

# Read all the vtk files. Ensure that they're read in chronological order.
vtk_files = [
    os.path.join(root, datafile)

    for root, folders, files in os.walk(BASEDIR)
    for datafile in files
    if datafile.lower().endswith('vtk')
]
readTimeInstant = lambda x: float(re.findall('\d+\.?\d*[eE]?[-+]?\d*', x)[0])
vtk_files.sort(key=readTimeInstant)

print('{} vtk files have been found'.format(len(vtk_files)))
print()


print('- - - - - Loading VTK surfaces - - - - -')

# Variables to be extracted
all_points = []
time = []
zmax_list = []
zmax_list2 = []

for i, filename in enumerate(vtk_files):
    print('Loading files... ', end='')
    print('{}/{}'.format(i+1, len(vtk_files)))

    # Read the file
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(filename)
    reader.Update()
    data = reader.GetOutput()

    # - Here is where the magic happens - #

    # Read the file:
    time.append(readTimeInstant(filename))

    # Point-By-Point analysis [slower]
    surface_points = np.array(
        [data.GetPoint(i) for i in range(data.GetNumberOfPoints())]
    )
        ## Extract the maximum index
    xind, yind, zind = 0, 1, 2  # <-- Use whatever index is down for your case
    zmax_list.append( surface_points.max(axis=0)[zind] )

    # Evaluate only bounds [faster]
    bbox = data.GetBounds()
    xmin, xmax, ymin, ymax, zmin, zmax = bbox
    zmax_list2.append(zmax)

print()

print('- - - - - Writting Results  - - - - -')
with open(output_file, 'w') as f:
    f.write('Time[s]:\tZ Max[m]:\n')
    for t,z in zip(time, zmax_list):
        f.write('{} \t{}\n'.format(t, z))
print('results have been saved into "{}"' \
      .format(os.path.abspath(output_file)))
print()

# -- Both graphs shall be equal --
print('- - - - - Plotting Results  - - - - -')
plt.plot(time, zmax_list, label='zmax_list')
plt.plot(time, zmax_list2, label='zmax_list2')
plt.xlabel('Time [s]:')
plt.ylabel('Height [m]:')
plt.grid()
plt.legend()
plt.show()
print()

print('DONE!')

I've tried the code with the sloshing3D tank and given me this graph (https://i.imgur.com/GR6rrHu.png), so it I hope it works for you as well.

Hopes it helps!


All times are GMT -4. The time now is 07:44.