動作環境

GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)

Numpy / Matplotlib > Meshgridでない座標をMeshgridに変換する > v0.1
の続き。

This article is related to ADDA (light scattering simulator based on the discrete dipole approximation).

ADDAの出力ファイル（dipoleのx,y,z座標)を表示してみた。

入力ファイル

ADDAの出力ファイル（dipoleのx,y,z座標) を含むIntField-Y。

$ head IntField-Y 
x y z |E|^2 Ex.r Ex.i Ey.r Ey.i Ez.r Ez.i
-0.2166615618 -1.516630933 -5.416539045 0.5516786267 -0.02931061187 0.02652841941 0.1694307002 0.7044162375 -0.09245370677 -0.1290700275
0.2166615618 -1.516630933 -5.416539045 0.5516786268 0.0293106112 -0.02652841834 0.1694306999 0.7044162377 -0.09245370673 -0.1290700275
-1.083307809 -1.083307809 -5.416539045 0.5009708155 -0.1210777782 0.1210667177 0.3201499555 0.5959123603 -0.06902233265 -0.09634427508
-0.6499846854 -1.083307809 -5.416539045 0.5264383492 -0.11184864 0.05898299433 0.1712903828 0.6818062435 -0.06862119958 -0.1074254128
-0.2166615618 -1.083307809 -5.416539045 0.7409394521 -0.04313270512 0.006539303629 0.1179657327 0.8426239137 -0.06420657934 -0.1047988574
0.2166615618 -1.083307809 -5.416539045 0.7409394524 0.04313270474 -0.00653930278 0.1179657325 0.8426239139 -0.06420657931 -0.1047988574
0.6499846854 -1.083307809 -5.416539045 0.5264383496 0.1118486397 -0.0589829935 0.1712903821 0.6818062441 -0.0686211995 -0.1074254129
1.083307809 -1.083307809 -5.416539045 0.5009708155 0.1210777777 -0.121066717 0.3201499543 0.5959123613 -0.06902233252 -0.09634427512
-1.083307809 -0.6499846854 -5.416539045 0.6583579786 -0.0961687486 0.04927804902 0.3522909973 0.7189310502 -0.04061949379 -0.0637218901

code

Jupyter code.

showChebyshev_170909.ipynb

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from pylab import rcParams
import sys


'''
v0.1 Sep. 09, 2017
   - read [IntField-Y] file

=== based on [toMeshgrid_170902.ipynb] ===

v0.1 Sep. 02, 2017
   - add get_meshgrid_from_xyzArray()
   - add func_z()
   - display 3D surface with lines
      + ref: https://stackoverflow.com/questions/30497737/applying-colormaps-to-custom-axis-in-matplotlib-3d-surface
'''

# coding rule: PEP8


rcParams['figure.figsize'] = 15, 10


def get_meshgrid_from_xyzArray(xar, yar, zar):
    # mx, my, mz : in meshgrid
    #
    xuniq = np.unique(xar)
    yuniq = np.unique(yar)
    mz = np.empty([len(yuniq), len(xuniq)])
    for ix in range(len(xuniq)):
        for iy in range(len(yuniq)):
            xx, yy = xuniq[ix], yuniq[iy]
            for idx in range(len(xar)):
                tx, ty = xar[idx], yar[idx]
                if abs(tx - xx) >= sys.float_info.epsilon:
                    continue
                if abs(ty - yy) >= sys.float_info.epsilon:
                    continue
                mz[iy][ix] = zar[idx]
    mx, my = np.meshgrid(xuniq, yuniq)
    return mx, my, mz

INPFILE = 'IntField-Y'

dat = np.genfromtxt(INPFILE, delimiter=' ', skip_header=1)

xpar, ypar, zpar = [], [], []
for elem in dat:
    axp, ayp, azp = elem[:3]
    xpar += [axp]
    ypar += [ayp]
    zpar += [azp]
    

# 2. from x,y,z arrays
res = get_meshgrid_from_xyzArray(xpar, ypar, zpar)
gx2, gy2, gz2 = res
ax2 = plt.subplot2grid((1, 1), (0, 0), projection='3d')
surf2 = ax2.plot_surface(gx2, gy2, gz2, shade=False,
                         facecolors=plt.cm.Set2((gx2-gx2.min())/(gx2.max()-gx2.min()))
                         )
plt.draw()

#   draw lines on the surface
lines = np.array(surf2.get_edgecolor())
# make lines white, and keep alpha==1. It's an array of colors like this: [r,g,b,alpha]
surf2.set_edgecolor(lines * np.array([0, 0, 0, 0]) + 1)

plt.show()

結果 > 球形粒子

粒子の形状を球形とした場合。

半球だけの表示となった。

結果 > Chebyshev粒子

Chebyshevの場合（パラメータは忘れた)。

失敗

「mx, myが同じ位置のmzを拾う」というスタンスなので、半球のみとなってしまっているようだ。

XとZの交換

XとZの値を交換してみた。

# 2. from x,y,z arrays
#res = get_meshgrid_from_xyzArray(xpar, ypar, zpar)
res = get_meshgrid_from_xyzArray(zpar, ypar, xpar)
gx2, gy2, gz2 = res
...

球形粒子の座標なのでX,Z交換は等価であるはず。

半球の形状になっていない。

前向きの考え
- 計算効率化のためにX,Y,Z座標のとり方が特殊になっている。
後向きの考え
- ソフトのバグ
  - 7of9作成部分だった。。。