Matplotlibによる応力・変形図の図化出力（3D対応）

前回紹介したMatplotlibによる図化出力ですが、以下の手順で3次元トラスの応力解析結果の図化出力にも拡張できました。蛇足ですが説明します。

図化出力のプログラムの頭に以下のコードを記載します。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

また、グラフのスクリプトを記述する際には、以下を記載します。

ax1 = fig.add_subplot(2,1,1, projection='3d')

以上を考慮するとMatplotlibによるグラフにおいてz軸の概念が導入され、2次元のグラフ描画と同様のルールでグラフの表示カスタマイズ（軸最大値の設定、表示する角度の設定など）も可能となります。Matplotlibによる3次元グラフ描画の詳細については以下のサイトに記載があります。

https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html

上記を踏まえて組んだプログラムにより、以下のように3次元のグラフが描画できました。意外といろいろなことができるんですね。Matplotlib。

参考として、今回の図化出力に用いた3次元トラスの応力解析結果の描画部分のソースコードを以下に示します。

#変位図、応力図を描く    
fig = plt.figure(figsize=(6, 9))
max_x=5
max_y=5#グラフの表示範囲を設定
max_z=5

ax1 = fig.add_subplot(2,1,1, projection='3d')
ax2 = fig.add_subplot(2,1,2, projection='3d')
for e in getEleTags():
    color_tag=0
    nodes = eleNodes(e)
    Nnodes = len(nodes)

    xyz = np.zeros((Nnodes,3), dtype=np.double)
    uu   = np.zeros((Nnodes,3), dtype=np.double)
    for i in range(Nnodes):
        xyz[i,:] = nodeCoord(nodes[i])#return coordinate
        uu[i,:] = nodeDisp(nodes[i])#return current displacement

    x = np.zeros(12,dtype=np.double)
    y = np.zeros(12,dtype=np.double)
    z = np.zeros(12,dtype=np.double)    
    u = np.zeros(12,dtype=np.double)
    v = np.zeros(12,dtype=np.double)
    t = np.zeros(12,dtype=np.double) 

    conec = [0, 1]

    x = xyz[conec,0]
    y = xyz[conec,1]
    z = xyz[conec,2]
    x[2::3] = np.nan
    y[2::3] = np.nan
    z[2::3] = np.nan

    mag_factor = 20.#変形図の拡大率
    u = xyz[conec,0] + mag_factor*uu[conec,0]
    v = xyz[conec,1] + mag_factor*uu[conec,1]
    t = xyz[conec,2] + mag_factor*uu[conec,2]
    u[2::3] = np.nan
    v[2::3] = np.nan
    t[2::3] = np.nan
    
    temp=basicForce(e)#各部材の応力をプロット
    if temp[0] >=0:
        if max_tens == 0:
            color_tag=temp[0]
        else:
            color_tag=temp[0]/max_tens*1       
        ax1.text(np.mean(x),np.mean(y),np.mean(z),round(temp[0],2))
        ax1.plot(x,y,z,color=[0,float(1-color_tag),1],lw=area[e-1])
    else:
        if max_comp == 0:
            color_tag=temp[0]
        else:        
            color_tag=-temp[0]/max_comp*1
        ax1.text(np.mean(x),np.mean(y),np.mean(z),round(temp[0],2))
        ax1.plot(x,y,z,color=[1,0,float(1-color_tag)],lw=area[e-1])

    ax2.plot(x,y,z,marker="o",color="0.5")       
    ax2.plot(u,v,t,color="0.",linestyle="dashed")#変形図をプロット

#各グラフのフォーマットの設定
ax1.set_title("Stress Diagram")
ax1.view_init(elev=40, azim=180)
ax1.set_xlim(-1,max_x)
ax1.set_ylim(-1,max_y)
ax1.set_zlim(-1,max_z)
ax1.grid(False)
ax1.set_xticks([]) 
ax1.set_yticks([]) 
ax1.set_zticks([]) 
ax2.set_title("Displacement Diagram")
ax2.view_init(elev=40, azim=180)
ax2.set_xlim(-1,max_x)
ax2.set_ylim(-1,max_y)
ax2.set_zlim(-1,max_z)
ax2.grid(False)
ax2.set_xticks([]) 
ax2.set_yticks([]) 
ax2.set_zticks([]) 
plt.show()