Parabolic Bowl

更新空いてしまいました。

Parabolic Bowl

先週末から平面二次元の計算を始めました。 手始めにparabolic bowlの図化をpythonで書いてみました。 気が向いたら回してみてください。DBに置いておきます。 さすがに遅すぎるので実計算には使えそうにないです。動画っぽくしたいなら、ImageMagickでgifにするほうが良いかもしれないです。

ソースコードはこんな感じです。ほぼグラフのところで、解析解の計算は数行です。やっぱりI/Oは手間ですね。

# -*- coding: UTF-8 -*-
# use python3.4
## @package
#
#
#
###%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

class MyGragh:
    def __init__(self,X,Y,zb,wL,idX,idY):
        self.X,self.Y,self.idX,self.idY = X,Y,idX,idY
    #make gragh
        self.fig = plt.figure(figsize=(12, 10))
    #make 3Dgragh
        self.ax = self.fig.add_subplot(2, 1, 1, projection='3d')
        self.ax.set_zlim3d(0, 14)
        self.ax.view_init(elev=33, azim=38)
        Xm, Ym = np.meshgrid(self.X, self.Y)
        h = list(map(list,zip(*zb)))
        self.ax.plot_surface(Xm, Ym, h, color='k', zorder=1, rstride=1, cstride=1,linewidth=0.1, antialiased=False,alpha=0.3)
        h = list(map(list,zip(*wL)))
        self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1)
    #make 2Dgragh
        self.ax2 = self.fig.add_subplot(2, 2, 3)
        plt.title('X = '+str(self.X[self.idX]))
        h = zb[self.idX,:]
        plt.plot(self.Y,h, color='k')
        h = wL[self.idX,:]
        self.lineX = self.ax2.plot(self.Y,h, color='b')

        self.ax3 = self.fig.add_subplot(2, 2, 4)
        plt.title('Y = '+str(self.Y[self.idY]))
        h = zb[:,self.idY]
        plt.plot(self.X,h, color='k')
        h = wL[:,self.idX]
        self.lineY = self.ax3.plot(self.X,h, color='b')

    def plotWL(self,wL,t,u,v):
    #make 3Dgragh
        self.surf.remove()
        Xm, Ym = np.meshgrid(self.X, self.Y)
        h = list(map(list,zip(*wL)))
        self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1)
    #make 2Dgragh
        for c in self.lineX : c.remove()
        h = wL[self.idX,]
        self.lineX = self.ax2.plot(self.Y,h, color='b')
        for c in self.lineY : c.remove()
        h = wL[:,self.idY]
        self.lineY = self.ax3.plot(self.X,h, color='b')

    #make title
        self.fig.suptitle('t='+str(t)+'[sec]'+' U='+str(round(u,3))+'[m/s]'+' V='+str(round(v,3))+'[m/s]')
        plt.pause(0.001)

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        plt.savefig('output.png')
        return False

# main
nx, ny = 21, 21
X = np.arange(-10, 11, 1)
Y = np.arange(-10, 11, 1)
a = 7
h0 = 4
zb = np.array([[h0/a**2*(X[i]**2+Y[j]**2) for j in range(ny)] for i in range(nx)])
eta = a/3
tPriod=60
omega = 2*np.pi/tPriod
#wL = lambda t: np.array([[ eta*h0/a**2*(2*X[i]*np.cos(t*omega)+2*Y[j]*np.sin(t*omega)-eta)+h0 for j in range(ny)] for i in range(nx)])
wL = lambda t, x, y: eta*h0/a**2*(2*x*np.cos(t*omega)+2*y*np.sin(t*omega)-eta)+h0
wLr = lambda t: np.array([[ wL(t,X[i],Y[j]) if wL(t,X[i],Y[j]) - zb[i,j] > 0 else zb[i,j] for j in range(ny)] for i in range(nx)])
u = lambda t: -eta*omega*np.sin(t*omega)
v = lambda t: -eta*omega*np.cos(t*omega)

idX=10
idY=10
tmax = 600
with MyGragh(X,Y,zb,wLr(0),idX,idY) as g:
    for t in range(1,tmax):
        if t%10==0 : g.plotWL(wLr(t),t,u(t),v(t))

アジョイント法

とりあえず、矩形断面ではできたっぽくて、テンションが上がってます。パッと見せれる絵がないので勉強会のときにでも話します。