カテナリー曲線は，自重に対しては曲げモーメントの発生しない幾何学曲線として知られており，サグラダ・ファミリアやゲートウェイ・アーチの構造にも取り入れられいている。ここでは，実際にカテナリー曲線に対して自重を作用させ，曲げモーメントが本当に0になるかどうかを考察するようなプログラムを作ってみよう。

カテナリー曲線は次式で定式化できることが知られている。

$$ y=a\cosh \left(\dfrac{x}{a}\right) $$

定数項$a$に加えて，アーチの分割数$n$およびスパン$S$を入力変数としてカテナリーアーチを構成する線分群を出力するGHPythonスクリプトを記述すると次のようになる。Rhinoビューポートは3次元空間なので，ここでは$y=0$

として$x-z$平面にカテナリーアーチを作成している。なお，$\cosh$を計算する関数は組み込み関数として存在しないため，import mathでmathライブラリを利用する。

import ghpythonlib.components as gc
import math

dx=S/float(n)
x=[]
for i in range(n+1):
    x.append(-S/2.0+dx*i)

z=[]

for i in range(n+1):
    z.append(a*math.cosh(x[i]/a))

lines=[]
for i in range(n):
    p1=gc.ConstructPoint(x[i],0,z[i])
    p2=gc.ConstructPoint(x[i+1],0,z[i+1])
    lines.append(gc.Line(p1,p2))

$n$の値が十分大きくなったとき，線分の集合はカテナリー曲線に近づく。

左：$n=4$　右：$n=50$

この線分をOpenSees for Grasshopperで読み込み，両端をピン支持として自重を載荷することを考えよう。数学の世界ではカテナリー曲線は下に凸の曲線として定義されているが，建築構造に用いる場合，符号を反転させて上に凸のアーチとして用いることも多いので，ParamsタブのValue Listを右クリックしてValue Cycleに切り替えたものを用いて$z$座標の正負を切り替える仕組みを作っておく。以下のコードでは，$+1$と$-1$のいずれかを取るpmという入力変数を付加し，$z$座標に乗じている。

これでボタン一つで関数の下に凸／上に凸を切り替えられるようになったが，このままでは，$S$や$a$の値を変化させると，頂点$(0,a)$を一定にしたままアーチの端部が動いてしまう。アーチの端部は常に$z$座標が0の位置にある方が構造解析モデルとしては扱いやすいので，下に凸／上に凸で場合分けをして，それぞれの$z$座標を集めたベクトル$\boldsymbol{z}$の最大値／最小値を$\boldsymbol{z}$から引き算することで，端部を常に$z=0$の位置へと調整する。

import ghpythonlib.components as gc
import math

dx=S/float(n)
x=[]
for i in range(n+1):
    x.append(-S/2.0+dx*i)

z=[]

for i in range(n+1):
    z.append(a*math.cosh(x[i]/a)*pm)

zmin,zmax=min(z),max(z)
if pm==1:
    for i in range(n+1):
        z[i]-=zmax
else:
    for i in range(n+1):
        z[i]-=zmin

lines=[]
for i in range(n):
    p1=gc.ConstructPoint(x[i],0,z[i])
    p2=gc.ConstructPoint(x[i+1],0,z[i+1])
    lines.append(gc.Line(p1,p2))

以上でジオメトリは完成したので，AssembpleGeometriesを介してVisualizeModelに接続してみる。すると，指定された分割数で線分が認識されていることが確認できる。

続いて境界条件を設定する。境界条件はBoundaryコンポーネントでも生成できるが，ここではGHPython内で直接生成することにする。境界条件の入力データは，

$$ \small\begin{bmatrix}\begin{bmatrix}\textrm{節点番号}&X\textrm{方向の拘束の有無(0 or 1)}&Y\textrm{方向の拘束の有無(0 or 1)}&Z\textrm{方向の拘束の有無(0 or 1)}&X\textrm{軸回りの回転の拘束の有無(0 or 1)}&Y\textrm{軸回りの回転の拘束の有無(0 or 1)}&Z\textrm{軸回りの回転の拘束の有無(0 or 1)}\end{bmatrix},\cdots\end{bmatrix} $$

という形式で，自由度を拘束する場合には1，そうでない場合には0を入れる。今回のアーチは，$0$番目の節点と$n$番目の節点を拘束すれば良いので，コードは次のようになる。なお，面外方向への倒れを防ぐため，ここでは$X，Y，Z$方向の移動の拘束に加えて$Y$軸回りの回転も拘束している。

import ghpythonlib.components as gc
import math

dx=S/float(n)
x=[]
for i in range(n+1):
    x.append(-S/2.0+dx*i)

z=[]

for i in range(n+1):
    z.append(a*math.cosh(x[i]/a)*pm)

zmin,zmax=min(z),max(z)
if pm==1:
    for i in range(n+1):
        z[i]-=zmax
else:
    for i in range(n+1):
        z[i]-=zmin

lines=[]
for i in range(n):
    p1=gc.ConstructPoint(x[i],0,z[i])
    p2=gc.ConstructPoint(x[i+1],0,z[i+1])
    lines.append(gc.Line(p1,p2))

bounds=[
[0,1,1,1,1,0,0],
[n,1,1,1,1,0,0]
]

なお，このコードではboundsに境界条件の情報を2重リストで定義しているが，GHPythonの2重リストをOpenSees for Grasshopperで読み込むためには，Treeに変換する必要があるため，SetsタブのListグループにあるLTConvを用いてGHPythonの2重リストからTreeに変換をした上で，VisualizeModelに接続しよう。ここまでの操作が正しく行われていれば，両端部に境界条件が生成されていることが確認できる。