Commit f409c654 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

draw velocity

parent e03733a6
from math import *
import matplotlib
from matplotlib.tri import Triangulation
from matplotlib import pyplot
from matplotlib.mlab import griddata
import matplotlib.cm as cm
import numpy
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
R = 3.3e-3
V = -0.00536667580007
a = 154e-6
phi = .2
#R = 3.7e-3
#V = -1.7e-3
#a = 154e-6
#phi = .04
dt = 5e-2
n = 10
N = 300
vm = []
vml = []
vmt = []
m = []
vf = []
for i in range(N):
r = []
x = []
v = []
b1 = False
b2 = False
b3 = False
b4 = False
b5 = False
b6 = False
k = i+1
with open("metzger4/MetzgerB/particles_%05d.vtp"% k) as f :
for l in f :
if b1 and l=='</DataArray>\n':
b2 = True
if b1 and not b2:
w = l.split()
for ww in w:
r.append(float(ww))
if l == '<DataArray Name="Radius" NumberOfComponents = "1" type="Float32" format="ascii">\n':
b1 = True
if b3 and l=='</DataArray>\n':
b4 = True
if b3 and not b4:
w = l.split()
for ww in w:
x.append(float(ww))
if l == '<DataArray NumberOfComponents="3" type="Float32" format="ascii">\n':
b3 = True
if b5 and l=='</DataArray>\n':
b6 = True
if b5 and not b6:
w = l.split()
for ww in w:
v.append(float(ww))
if l == '<DataArray Name="Velocity" NumberOfComponents = "3" type="Float32" format="ascii">\n':
b5 = True
r = numpy.array(r)
x = numpy.array(x)
v = numpy.array(v)
dl = []
vl = 0
vl1 = 0
xl = 0
yl = 0
zl = 0
ii = 0
mmm = 10
for j in range(len(r)):
xl += x[3*j+0]
yl += x[3*j+1]
if x[3*j+1]<mmm:
mmm = x[3*j+1]
zl += x[3*j+2]
vf.append(mmm)
xl /= (len(v)/3.)
yl /= len(r)
zl /= (len(v)/3.)
hl = 0
for j in range(len(r)):
dl.append(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5)
if x[3*j+1]<yl:
hl = max(hl,yl-x[3*j+1])
dl1 = 0
dl2 = 0
dl3 = 0
dl4 = 0
for j in range(len(r)):
dl.append(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5)
if x[3*j]<0 and x[3*j+2]<0:
dl1 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl1)
if x[3*j]>0 and x[3*j+2]<0:
dl2 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl2)
if x[3*j]<0 and x[3*j+2]>0:
dl3 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl3)
if x[3*j]>0 and x[3*j+2]>0:
dl4 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl4)
dl.sort()
dl = numpy.array(dl)
# m.append((dl1+dl2+dl3+dl4)/4)
m.append(dl[int(.8*len(r))])
yy = mmm+m[0]
for j in range(len(r)):
if ((xl-x[3*j+0])**2+(yy-x[3*j+1])**2+(zl-x[3*j+2])**2)**.5<3.3e-3:
ii += 1
vl += v[3*j+1]
else:
vl1 += v[3*j+1]
vl /= ii
vl1 /= (len(r)-ii+1)
vm.append(vl)
vml.append(vl1/vm[0])
vmt.append((ii*vl/vm[0]+(len(r)-ii)*vl1/vm[0])/len(r))
print(ii)
vf = numpy.array(vf)
for i in range(len(vf)-1):
vf[i] = (vf[i+1]-vf[i])/dt
vf = vf[:-1]/vf[0]
tf = numpy.arange(0., N-1, 1)*dt*n/(m[0]/(-vf[0]))
t = numpy.arange(0., N, 1)*dt*n/(m[0]/(-vm[0]))
print(m[0])
############2
vm1 = []
m1 = []
for i in range(N):
r = []
x = []
v = []
b1 = False
b2 = False
b3 = False
b4 = False
b5 = False
b6 = False
k = i+1
with open("metzger5/MetzgerB/particles_%05d.vtp"% k) as f :
for l in f :
if b1 and l=='</DataArray>\n':
b2 = True
if b1 and not b2:
w = l.split()
for ww in w:
r.append(float(ww))
if l == '<DataArray Name="Radius" NumberOfComponents = "1" type="Float32" format="ascii">\n':
b1 = True
if b3 and l=='</DataArray>\n':
b4 = True
if b3 and not b4:
w = l.split()
for ww in w:
x.append(float(ww))
if l == '<DataArray NumberOfComponents="3" type="Float32" format="ascii">\n':
b3 = True
if b5 and l=='</DataArray>\n':
b6 = True
if b5 and not b6:
w = l.split()
for ww in w:
v.append(float(ww))
if l == '<DataArray Name="Velocity" NumberOfComponents = "3" type="Float32" format="ascii">\n':
b5 = True
r = numpy.array(r)
x = numpy.array(x)
v = numpy.array(v)
dl = []
vl = 0
vl1 = 0
xl = 0
yl = 0
zl = 0
ii = 0
mmm = 10
for j in range(len(r)):
xl += x[3*j+0]
yl += x[3*j+1]
if x[3*j+1]<mmm:
mmm = x[3*j+1]
zl += x[3*j+2]
xl /= (len(v)/3.)
yl /= len(r)
zl /= (len(v)/3.)
hl = 0
for j in range(len(r)):
dl.append(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5)
if x[3*j+1]<yl:
hl = max(hl,yl-x[3*j+1])
dl1 = 0
dl2 = 0
dl3 = 0
dl4 = 0
for j in range(len(r)):
dl.append(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5)
if x[3*j]<0 and x[3*j+2]<0:
dl1 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl1)
if x[3*j]>0 and x[3*j+2]<0:
dl2 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl2)
if x[3*j]<0 and x[3*j+2]>0:
dl3 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl3)
if x[3*j]>0 and x[3*j+2]>0:
dl4 = max(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5,dl4)
dl.sort()
dl = numpy.array(dl)
# m1.append((dl1+dl2+dl3+dl4)/4)
m1.append(dl[int(.8*len(r))])
yy = mmm+m1[0]
for j in range(len(r)):
if ((xl-x[3*j+0])**2+(yy-x[3*j+1])**2+(zl-x[3*j+2])**2)**.5<3.3e-3:
ii += 1
vl += v[3*j+1]
else:
vl1 += v[3*j+1]
vl /= ii
vl1 /= (len(r)-ii+1)
vm1.append(vl)
vml.append(vl1/vm1[0])
vmt.append((ii*vl/vm1[0]+(len(r)-ii)*vl1/vm1[0])/len(r))
print(ii)
t1 = numpy.arange(0., N, 1)*dt*n/(m1[0]/(-vm1[0]))
print(vm1[0])
vm0 = (4*numpy.array(vm)-numpy.array(vm1))/3.
m0 = (4*numpy.array(m)-numpy.array(m1))/3.
t0 = numpy.arange(0., N, 1)*dt*n/(m0[0]/(-vm0[0]))
vm1 = numpy.array(vm1)/(vm1[0])
m1 = numpy.array(m1)/m1[0]
vm = numpy.array(vm)/(vm[0])
m = numpy.array(m)/m[0]
vm0 = numpy.array(vm0)/(vm0[0])
m0 = numpy.array(m0)/m0[0]
#####################3
#vm2 = []
#m2 = []
#for i in range(N):
# r = []
# x = []
# v = []
# b1 = False
# b2 = False
# b3 = False
# b4 = False
# b5 = False
# b6 = False
# k = i+1
# with open("metzger6/MetzgerB/particles_%05d.vtp"% k) as f :
# for l in f :
# if b1 and l=='</DataArray>\n':
# b2 = True
#
# if b1 and not b2:
# w = l.split()
# for ww in w:
# r.append(float(ww))
#
# if l == '<DataArray Name="Radius" NumberOfComponents = "1" type="Float32" format="ascii">\n':
# b1 = True
#
#
#
# if b3 and l=='</DataArray>\n':
# b4 = True
#
# if b3 and not b4:
# w = l.split()
# for ww in w:
# x.append(float(ww))
#
# if l == '<DataArray NumberOfComponents="3" type="Float32" format="ascii">\n':
# b3 = True
#
#
#
# if b5 and l=='</DataArray>\n':
# b6 = True
#
# if b5 and not b6:
# w = l.split()
# for ww in w:
# v.append(float(ww))
#
# if l == '<DataArray Name="Velocity" NumberOfComponents = "3" type="Float32" format="ascii">\n':
# b5 = True
#
#
# r = numpy.array(r)
# x = numpy.array(x)
# v = numpy.array(v)
# dl = []
# vl = 0
# vl1 = 0
# xl = 0
# yl = 0
# zl = 0
# ii = 0
# mmm = 10
# for j in range(len(r)):
# xl += x[3*j+0]
# yl += x[3*j+1]
# if x[3*j+1]<mmm:
# mmm = x[3*j+1]
# zl += x[3*j+2]
#
# xl /= (len(v)/3.)
# yl /= len(r)
# zl /= (len(v)/3.)
# hl = 0
# for j in range(len(r)):
# dl.append(((xl-x[3*j+0])**2+(zl-x[3*j+2])**2)**.5)
# if x[3*j+1]<yl:
# hl = max(hl,yl-x[3*j+1])
#
# dl.sort()
# dl = numpy.array(dl)
#
# m2.append(dl[int(0.95*len(r))])
#
# RR = mmm+2*dl[int(.95*len(r))]
# for j in range(len(r)):
# if x[3*j+1]<RR:
# ii += 1
# vl += v[3*j+1]
# else:
# vl1 += v[3*j+1]
# vl /= ii
# vl1 /= (len(r)-ii)
#
# vm2.append(vl)
# print(vl1/vm2[0],vl/vm2[0],(ii*vl/vm2[0]+(len(r)-ii)*vl1/vm2[0])/len(r))
#t2 = numpy.arange(0., N, 1)*dt*n/(m2[0]/(-vm2[0]))
#print(vm2[0])
#vm2 = numpy.array(vm2)/(vm2[0])
#m2 = numpy.array(m2)/m2[0]
fig, axarr = pyplot.subplots(2,sharex=True)
axarr[0].plot(t,m,label='$\epsilon/2$')
axarr[0].plot(t1,m1,label='$\epsilon$')
axarr[0].plot(t0,m0,'k--',label='$\epsilon^*$')
legend = axarr[0].legend(loc='upper center', bbox_to_anchor=(0.5, 1.25),fancybox=True, shadow=True, ncol=3)
frame = legend.get_frame()
frame.set_facecolor('0.90')
for label in legend.get_texts():
label.set_fontsize('large')
for label in legend.get_lines():
label.set_linewidth(1.2)
axarr[0].set_xlim([0,250])
axarr[0].set_ylabel('$R^*$')
axarr[0].grid(True)
#axarr[1].plot(tf,vf,label='$\epsilon=1.8e-3$')
#axarr[1].plot(t,vm1,label='$\epsilon=1.8e-3$')
#axarr[1].plot(t0,vm0,'b--',label='$\epsilon=1.8e-3$')
axarr[1].plot(tf,vf,label='$\epsilon=1.8e-3$')
axarr[1].set_xlim([0,250])
#axarr[0].set_ylim([.9,1.5])
axarr[1].set_xlabel('$T^*$')
axarr[1].set_ylabel('$V^*$')
axarr[1].grid(True)
#pyplot.show()
fig.savefig('MetzgerBeps.png',bbox_inches='tight',dpi=800)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment