【FreeCAD使い方講座】FreeCADのmatplotlibを使ってGyroid関数を3次元プロットしてみる

以前の記事で、

を簡単に使ってみる方法を紹介したことがあります。
これにより、FreeCADユーザーはわざわざmatlibplot専用のPython環境を整えてあげなくても、FreeCADを立ち上げる感覚そのままにnumpyやscipyから提供される高度な関数をある程度利用することが可能です。

【FreeCAD初心者ガイド】FreeCAD付属のmatplotlibを使って3次元モデリングに便利なプロトタイプ関数を思索する

FreeCADに標準で組込まれているPython製高機能プロッター・matplotlibの基本的な使い方を紹介

今回の記事では、matplotlibを使ったPythonコードの活用の一例として、以下の記事で紹介していた

「Gyroid関数」

の解析解を導いて3次プロットしてみるまでを紹介します。

【pythonで動くG-codeビルダー】G-coordinator/gcoordinatorをLinuxにインストール＆動作確認してみる

Pythonソースコードで3Dプリンター向けのgcodeが生成できる・『G-coordinator』をLinuxへ導入します。

Gyroid関数の解析解

まずGyroid関数は以下のような陰関数で与えられます。

\displaystyle { sin(x) cos(y) + sin(y) cos(z) + sin(z) cos(x) = 0 }

Eq. (gyroid-1)

Gyroidは陰関数で書くとサイクリックな項の和として表せるため、綺麗な見た目をしています。
ただし、このまま３次元グラフを描く際に、Mathematicaのような高度な数学的な処理が可能なソフトウェアが必要になります。
今回はmatplotlibから3次元プロットを試みたいのですが、数値解的な手法を使うと、ソルバをどう設定すれば精度が良いか...など他に考慮しなければならない点が多くなります。
こういったことから、matplotlibで３次元グラフ化するためには、

解析解(陽関数)

を得ることが出来ないかを最初に考えると良いでしょう。
式変形した結果だけを示しますが、Gyroid関数は以下のような2変数関数として解析解を持ちます。

\displaystyle { z = \arcsin \biggl\{ \frac{-\cos{x}\sin{y}\cos{y} \pm \sin{x}\sqrt{\sin^2{x} + \cos^2{y} - \cos^2{x} \sin^2{y}}}{\sin^2{x} + \cos^2{y}} \biggr\} }

Eq. (gyroid-2)

陰関数形とは違い割と綺麗ではないですが、計算としてはより単純です。
注意が必要なのは、分子の中にある平方根の取り扱いで、ルートの中身が負の値をとるとき、$$z$$は実数解を持ちません。
また、二次方程式由来の$$\pm$$からも分かる通り、Gyroid関数の解析解は「2枚張り」の対で存在しています。
matplotlibで3次元描画する際にはこの点を注意してスクリプトを書いていきましょう。

matplotlibでの描画

では先程の解析解の考え方を基にして、Pythonコードに落とし込んでいきます。

開く/閉じる

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

#--- メッシュグリッドの生成(1回目) ---

N = 1000
R1 = -np.pi/2
R2 = np.pi/2
x = np.linspace(R1, R2, N)
y = np.linspace(R1, R2, N)
X, Y = np.meshgrid(x, y)

#--- メッシュグリッドの生成(2回目) ---

A1 = np.sin(X)**2 + np.cos(Y)**2
A2 = np.cos(X)*np.sin(Y)
Z0 = A1 - A2**2
X,Y = np.where(Z0 >= 0, X, np.nan), np.where(Z0 >= 0, Y, np.nan)

#--- Gyroid関数の計算 ---

A1 = np.sin(X)**2 + np.cos(Y)**2
A2 = np.cos(X)*np.sin(Y)
A3 = (-A2*np.cos(Y) - np.sin(X)*np.sqrt(A1 - A2**2))/A1
A4 = (-A2*np.cos(Y) + np.sin(X)*np.sqrt(A1 - A2**2))/A1
Z1 = np.arcsin(A3)
Z2 = np.arcsin(A4)

#--- ワイヤーフレーム図として描画 ---

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X,Y,Z1)
ax.plot_wireframe(X,Y,Z2)

plt.show()

このスクリプトのポイントだけ掻い摘むと、matplotlibのワイヤーフレーム描画には、局面表示させるためのメッシュグリッドが必要になりますが、先程の節でも述べたように、実数解を持たない領域は避ける必要があります。
ということで、2回目のメッシュグリッド生成のときに、

np.where

を使って、有効なワイヤーフレームに矯正しています。
なおここでは、matplotlibではデータ点に

nan

がある場合には描画をスキップしてくれるを期待して、

np.nan

をワイヤーフレーム座標に入れています。
これを実行すると、
900x711

という感じのGyroid局面が描けます。
メッシュに裂け目があるのは、前述した通り、解析解が「2枚張り」であるがゆえの2つのワイヤーフレーム面が表示として見えています。
2つのワイヤーフレームの表示を上手く接続できる方法もあるかもしれませんが、そちらはmatplotlibのテクニックの問題になるため、これ以上は深い入りするのはやめておきます。

まとめ

以上、今回はGyroid関数の解析解を解いてみようというそれだけの企画でした。
というのも、以前やった「G-coordinator」の

サンプルコード

にあるGyroid関数の解き方が何かとトリッキーであると感じたためで、次回以降では今回の解析式を利用して座標を求めていくようにスクリプトの補正をしてみたいと思います。