CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

Postprocessing .vtk file OpenFOAM

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By Antimony

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 21, 2020, 12:06
Lightbulb Postprocessing .vtk file OpenFOAM
  #1
New Member
 
borA
Join Date: Dec 2019
Posts: 6
Rep Power: 2
kaboka is on a distinguished road
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;
   }  
);
kaboka is offline   Reply With Quote

Old   January 22, 2020, 01:29
Default
  #2
Senior Member
 
Join Date: Aug 2013
Posts: 388
Rep Power: 10
Antimony is on a distinguished road
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 likes this.
Antimony is offline   Reply With Quote

Old   January 22, 2020, 10:11
Default
  #3
New Member
 
borA
Join Date: Dec 2019
Posts: 6
Rep Power: 2
kaboka is on a distinguished road
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 is offline   Reply With Quote

Old   March 24, 2020, 04:56
Default
  #4
New Member
 
borA
Join Date: Dec 2019
Posts: 6
Rep Power: 2
kaboka is on a distinguished road
Quote:
Originally Posted by Antimony View Post
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?
kaboka is offline   Reply With Quote

Old   March 31, 2020, 23:18
Default
  #5
Senior Member
 
Join Date: Aug 2013
Posts: 388
Rep Power: 10
Antimony is on a distinguished road
Hi,

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

Cheers,
Antimony
Antimony is offline   Reply With Quote

Old   April 1, 2020, 07:05
Default
  #6
New Member
 
borA
Join Date: Dec 2019
Posts: 6
Rep Power: 2
kaboka is on a distinguished road
Quote:
Originally Posted by Antimony View Post
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
kaboka is offline   Reply With Quote

Old   April 1, 2020, 09:11
Default
  #7
Senior Member
 
zhangyan's Avatar
 
Yan Zhang
Join Date: May 2014
Posts: 106
Rep Power: 8
zhangyan is on a distinguished road
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 View Post
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
__________________
https://openfoam.top/en
zhangyan is offline   Reply With Quote

Old   May 21, 2020, 09:45
Post Possible Code solution?
  #8
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 34
Rep Power: 4
crubio.abujas is on a distinguished road
Hi kaboka,

Quote:
Originally Posted by kaboka View Post
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 View Post
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!
crubio.abujas is online now   Reply With Quote

Reply

Tags
freesurface, vtk

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[swak4Foam] groovyBC in openFOAM-2.0 for parabolic velocity bc ofslcm OpenFOAM Community Contributions 25 March 6, 2017 10:03
how to calculate mass flow rate on patches and summation of that during the run? immortality OpenFOAM Post-Processing 101 February 2, 2017 23:36
[swak4Foam] funkyDoCalc with OF2.3 massflow NiFl OpenFOAM Community Contributions 11 November 1, 2016 06:43
[OpenFOAM.org] Compile OF 2.3 on Mac OS X .... the patch gschaider OpenFOAM Installation 225 August 25, 2015 19:43
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51


All times are GMT -4. The time now is 08:21.