[kaboka]

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

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]

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]

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.

#!/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]

Originally Posted by Antimony

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]

Hi,

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

Cheers,
Antimony

[kaboka]

Originally Posted by Antimony

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:

OSError: coords not found. 


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

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:

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]

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

Originally Posted by Antimony

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]

Hi kaboka,

Originally Posted by kaboka

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:

OSError: coords 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:

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.

Originally Posted by kaboka

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

This is the website:

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.

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

[kaboka]

Hello Carlos,

Thanks a lot, that was super explained. I will work on it.

Bests

[mkuzaay]

for OpenFOAM-v2012+

freeSurface
     {
        type            surfaces;
        libs            (sampling);
        interpolationScheme cell;
        writeControl    runTime;
        writeInterval   0.1;
        surfaceFormat   vtk;
        fields
        (
            alpha.water
        );
        surfaces
        (
            freeSurface
            {
                type        isoSurfaceCell;
                isoField    alpha.gas;
                isoValue    0.5;
                interpolate true;
                regularise  false;
            }    
         );
     }