[kaboka]

Hello everyone,

I have sloshing case inside the 3D container. I got the values with .vtk code as follows. I wrote that code inside the controlDict and it worked.

PHP 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;
   }  
); 


Now based on this website:

https://openfoamwiki.net/index.php/Tip_Surface_elevation_in_time

I would like to plot the specific point in the free surface.

PHP 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")

# 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() 


I'm using this code but I am getting an OSError: coords not found.
Could anyone help me how to create the coordinate file and plot it in the python, please ?

Best regards