bin/plotForces.py

   1 #!/usr/bin/python
   2
   3 forcesfilename = 'forces/0/forces.dat'
   4
   5 import sys
   6 if len(sys.argv) != 1:
   7         print 'script assumes forces file ', forcesfilename
   8         sys.exit()
   9
  10 import re
  11 forceRegex=r"([0-9.Ee\-+]+)\s+\(+([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\)\s\(([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9 .Ee\-+]+)\)+\s\(+([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\)\s\(([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\s([0-9.Ee\-+]+)\)+"
  12 t = []
  13 fpx = []; fpy = []; fpz = []
  14 fvx = []; fvy = []; fvz = []
  15 mpx = []; mpy = []; mpz = []
  16 mvx = []; mvy = []; mvz = []
  17
  18 pipefile=open(forcesfilename,'r')
  19 lines = pipefile.readlines()
  20
  21 for line in lines:
  22         match=re.search(forceRegex,line)
  23         if match:
  24                 t.append(float(match.group(1)))
  25                 fpx.append(float(match.group(2)))
  26                 fpy.append(float(match.group(3)))
  27                 fpz.append(float(match.group(4)))
  28                 fvx.append(float(match.group(5)))
  29                 fvy.append(float(match.group(6)))
  30                 fvz.append(float(match.group(7)))
  31                 mpx.append(float(match.group(8)))
  32                 mpy.append(float(match.group(9)))
  33                 mpz.append(float(match.group(10)))
  34                 mvx.append(float(match.group(11)))
  35                 mvy.append(float(match.group(12)))
  36                 mvz.append(float(match.group(13)))
  37
  38
  39 # combine pressure and viscous forces and moments
  40 fx = fpx
  41 fy = fpy
  42 fz = fpz
  43
  44 mx = mpx
  45 my = mpy
  46 mz = mpz
  47
  48 for i in range(1,len(t)):
  49         fx[i] += fvx[i]
  50         fy[i] += fvy[i]
  51         fz[i] += fvz[i]
  52         mx[i] += mvx[i]
  53         my[i] += mvy[i]
  54         mz[i] += mvz[i]
  55
  56 # plot forces
  57 import pylab
  58 pylab.xlabel('iteration')
  59 pylab.ylabel('force (N)')
  60 pylab.grid(True)
  61 #
  62 pylab.plot(t,fx,'-',label="fx")
  63 pylab.plot(t,fy,'-',label="fy")
  64 pylab.plot(t,fz,'-',label="fz")
  65
  66 #pylab.plot(t,mx,'-',label="mx")
  67 #pylab.plot(t,my,'-',label="my")
  68 #pylab.plot(t,mz,'-',label="mz")
  69
  70 pylab.legend()
  71 pylab.show()