20191014

three.jsのpdb_loaderを使って、カーボンナノチューブを表示した。

やったこと

three.jsを使って、pdbファイルのカーボンナノチューブを表示した。
カーボンナノチューブのpdbファイルは、vmdというソフトを使って書き出した。

参考図書に、「はじめてのThree.js」を使ったけれど、pdb_loaderがうまく動かなかった。
githubのサンプル通りに、書いたのに(もちろん一部改変しなければいThree.jsのアドレスとかの所は改変した)表示されなかった。

どうも、この部分、https://github.com/josdirksen/learning-threejs/blob/85f85aeb137f0b01ba96c8b1bfd7019134e596ef/chapter-08/12-load-pdb.html#L88-L111
が動いてくれていない。
コードの中で、geometry.vertices.forEachとか書いてあるけど、「geometry以下にverticesがない」様な事を示唆するエラーが出てくる。(undefinedがどうとか)

ググったら
(https://threejs.org/docs/#examples/en/loaders/PDBLoader)
が関係ありそうだったので見た。

公式refereceではコールバック関数が一つ。
対する「はじめてのThree.js」では二つ。

書籍だから、たぶん情報が古い。昔はこれで動いたんだろう。

コールバック関数を一つにして再度チャレンジ。
この関数から、

  • コールバック関数名.geometryAtoms
  • コールバック関数名.geometryBonds
  • コールバック関数名.json

という感じの3つのデータが取り出せる。
しかし、jsonに原子と結合手のデータが全部入っていたので、そっちを使った。


変更前

  • 原子を元素名ごとに色分けしてsphereで表示してた
  • 結合手をsylinderで表示してた

変更後

  • 原子はすべて単色、
  • 結合手は、lineで書く。

L88-L111を、下記のように書き換えたら大体うまく動いた。
動くという意味では改善だが、原子の大きさや色などを反映しなくなった点で改悪。
var loader = new THREE.PDBLoader();
var group = new THREE.Object3D();
var json=loader.load("http://adoresu-wo-ireru.sakura.ne.jp/sblo_files/anatano-no-saito-no-namae/tekitouna_foruda_demo_tukure/nanotube.pdb", function (geometry) {
var json = geometry.json;
for (var j = 0; j < json.atoms.length; j += 1) {
var sphere = new THREE.SphereGeometry(0.1,8,8);
var material = new THREE.MeshBasicMaterial( {color: 0xffff00} );
var mesh = new THREE.Mesh(sphere, material);
mesh.position.set(json.atoms[j][0], json.atoms[j][1], json.atoms[j][2]);
group.add(mesh);
};
for (var j = 0; j < json.bonds.length; j += 1) {
var material2 = new THREE.MeshBasicMaterial( {color: 0x0000ff} );
eval("var ln" + j + "= new THREE.Geometry();");
// console.log(json.bonds.length)
ii=0;
json.bonds[j].forEach(function(element) {
if (ii != 2) {
var aaa = element
eval("ln" + j + ".vertices.push(new THREE.Vector3(json.atoms[aaa][0], json.atoms[aaa][1], json.atoms[aaa][2]));");
// console.log(aaa);
ii+=1;
};
var material4 = new THREE.LineBasicMaterial({color: 0x0000ff});
eval("var line" + j + "= new THREE.Line(ln" + j + ", material4 );");
eval("group.add( line" + j + " );");
});
};
scene.add(group);
});

結果はこちら

pdb_loader_test.png

上の画面の様なCNTがクルクル回るのが、自分のインターネットブラウザの環境からは見える。

ネット上で3Dが動かせている。その内インタラクティブにもできるようになるだろう。

最終目標は、あくまでも「web、地震シミュ、VR」インタラクティブは、視点の移動位でいいかも。

posted by yuchan at 08:00 | Comment(1) | javascript

20191007

断面二次モーメントを、座標点の配列から計算する(2)断面二次モーメント編

断面二次モーメント

やろうとしている事

断面二次モーメントを数値計算的に算出しようとしている。(先々月の続き。)

計算の仕方は二つ

  1. 長さ寸法と、公式に当てはめて計算をする。

     公式の公式集。書籍よりDVD版の方が明らかにお得だが高い。

    • メリット
      • 正確な値
    • ディメリット
      • 形状によっては式が用意されていない。
  2. 座標値を使い、数値計算的に求める

  3. メリット
    • 断面形状にかかわらず、座標値さえあれば計算可能。
  4. ディメリット
    • 点の間隔、カーブ、数値積分の手法などで誤差が生じる。
    • 座標値を準備しなければならない。

1の方法が普通のやり方。できるなら1をやるべき。

断面二次モーメントの公式

$$\int _{\rm A}y^2\rm\ dA$$

上は、断面二次モーメントの式。

下は、8×8の正方形の断面二次モーメントを計算する際の積分のイメージ。

20190918130448_2.png

数値計算で断面二次モーメントを求めるための、計算手順

断面二次モーメントまでの道のり
  1. 座標値を手に入れる

  2.  - 断面一次モーメントを計算する。

    • 断面一次モーメントから図心を計算する。
    • 図心が原点となるように座標値を整える。(省略)
  3. 断面二次モーメントを計算する。

1. 座標値を手に入れる 四角形のx,y座標を得る関数
def coorsquare(r=4, center=[0,0], n=10, closed=True):
rx1=-r/2.0+center[0]
rx2=r/2.0+center[0]
ry1=-r/2.0+center[1]
ry2=r/2.0+center[1]

x = [rx1 + i*(rx2-rx1)/float(n) for i in range(n)]
y = [ry2 for i in range(n)]
#1
x = x + [rx2 for i in range(n)]
y = y + [ry2 - i*(ry2-ry1)/float(n) for i in range(n)]
#2
x = x + [rx2 - i*(rx2-rx1)/float(n) for i in range(n)]
y = y + [ry1 for i in range(n)]
#3
x = x + [rx1 for i in range(n)]
y = y + [ry1 + i*(ry2-ry1)/float(n) for i in range(n)]
#4
if closed:
x = x + [rx1]
y = y + [ry2]

return x, y


x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
上の関数を使って作ったx,yを作図。
import matplotlib.pyplot as plt

fig,ax1 = plt.subplots(1)
ax1.plot(x,y)
ax1.scatter(x,y)

for i in range(len(x)):
ax1.text(x[i],y[i],"%s"%(i))

ax1.set_xlim(-4,4)
ax1.set_ylim(-4,4)
ax1.set_aspect("equal")

plt.show()

shikaku.png

2. 面積を計算→断面一次モーメント→図心を計算。
def area(x, y):
A = sum([(y[i]+y[i+1])/2.0*(x[i+1]-x[i]) for i in range(len(x)-1)])
return A

x,y = coorsquare(r=5, center=[0,0], n=20, closed=True)
print(area(x,y))

def FirstInrt(x, y):
FIy = sum([(y[i]+y[i+1])**2/8.0*(x[i+1]-x[i]) for i in range(len(x)-1)])
FIx = sum([(x[i]+x[i+1])**2/8.0*(y[i]-y[i+1]) for i in range(len(x)-1)])
return FIx, FIy

print(FirstInrt(x, y))


def CntCS(x, y):
A = area(x, y)
FIx, FIy = FirstInrt(x, y)
Cntx, Cnty = FIx/A, FIy/A
return Cntx, Cnty

print("center(%s,%s)"%(CntCS(x, y)))
3. 断面二次モーメントの計算(シンプソン則)

二次の式なので、台形則でなくシンプソン則で計算した。この関数が、この文章の一番重要な所

def simpson(x,y):#sekibun_1/3simpson(3point)
Cntx, Cnty = CntCS(x, y)
x = np.array(x)-Cntx
y = np.array(y)-Cnty
SecondInrty,SecondInrtx = 0,0

x1 = np.array(x)**3/3.0
y1 = np.array(y)**3/3.0

for i in range(len(x))[::2][:-1]:
a1 = float(x[i+1]-x[i])
# if は、0割対策。
if a1 !=0.0:
a1 = y1[i] * (x[i+2]-x[i])*(-x[i+2]-2.0*x[i]+3.0*x[i+1])/a1/6.0
a2 = float((x[i+1]-x[i])*(x[i+2]-x[i+1]))
if a2 !=0.0:
a2 = y1[i+1] * (x[i+2]-x[i])**3/a2/6.0
a3 = float(x[i+2]-x[i+1])
if a3 !=0.0:
a3 = y1[i+2] * (x[i+2]-x[i])*(2.0*x[i+2]+x[i]-3.0*x[i+1])/a3/6.0

SecondInrty = SecondInrty + a1 + a2 + a3

b1 = y[i+1]-y[i]
if b1 !=0:
b1 = x1[i]*(y[i+2]-y[i])*(-y[i+2]-2*y[i]+3*y[i+1])/b1/6.0
b2 = (y[i+1]-y[i])*(y[i+2]-y[i+1])
if b2 !=0:
b2 = x1[i+1]*(y[i+2]-y[i])**3/6.0/b2
b3 = y[i+2]-y[i+1]
if b3 !=0:
b3 = x1[i+2]*(y[i+2]-y[i])*(2*y[i+2]+y[i]-3*y[i+1])/b3/6.0

SecondInrtx = SecondInrtx - (b1 + b2 + b3)

return SecondInrtx, SecondInrty


# 断面が中心にない場合でも、図心を原点に持って行った時の値になるようにした。
x,y = coorsquare(r=5, center=[7,7], n=8, closed=True)
print("1/3simpson(3point)%s,%s"(simpson(x,y)))

# nは、偶数にしないと、正しく計算されない。
#うまくいかない。
x,y = coorsquare(r=5, center=[0,0], n=7, closed=True)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

# 全体で偶数になっていても、一辺が奇数ではだめ。
#うまくいかない
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)
y.insert(-1,(2.5+1.875)/2.0)
x.insert(-1,-2.5)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

# nが偶数なら、不等間隔でも計算される。
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)
x.insert(-9,-2.4)
y.insert(-9,-2.5)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

for ループは、2個飛ばしで行う。

3点使って三個飛ばしだと、間が空く。

3点使って二個飛ばしだと、間が埋まる。

ごめん、うまく書けない。こういう事。

figure_1.png

ってなっちゃうから、3個飛ばしにすればいいわけじゃない。

3個使うけど、最後の一個は次のループでも使うから、二個飛ばしが正解。
(ちなみに、pythonで二個飛ばしの書き方は、[::2]。)
こういうindexingの所で、いつも試行錯誤してる。頭だけで出来ないとだめだろうか?

figure_1-2.png

前のループでi+2だったところを、次のループではiとして扱っている。

積分は不等間隔のシンプソン則で模擬。不等間隔のシンプソン則は、ココを参考にして作った。

hutou_kankaku.png

一辺のサンプルが奇数ならば正しく計算できる。

ひし形(45度回転。一部不等間隔。)
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)x.insert(-9,-2.4)
y.insert(-9,-2.5)
xy= np.dot(np.array([x,y]).T, [[np.cos(1/4.0*np.pi), np.sin(1/4.0*np.pi)],[-np.sin(1/4.0*np.pi), np.cos(1/4.0*np.pi)]]).T
x=xy[0]
y=xy[1]

print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

hishi_gata.png

前回は、これが正しく計算できなかったが、今回は不等間隔シンプソン則のおかげで計算できた。
(サンプル数の偶奇に気を付けなければならないが)

穴あき(正方形)
x1,y1 = coorsquare(r=5, center=[0,0], n=8, closed=True)
x2,y2 = coorsquare(r=3, center=[0,0], n=4, closed=True)
x2=x2[::-1]
y2=y2[::-1]
Ix1,Iy1 = simpson(x1,y1)
Ix2,Iy2 = simpson(x2,y2)
a = Ix1+Ix2
b = Iy1+Iy2
print("1/3simpson(3point)%s,%s"%(a,b))
import matplotlib.pyplot as plt

fig,ax1 = plt.subplots(1)
ax1.plot(x1,y1)
ax1.plot(x2,y2)
ax1.scatter(x1,y1)
ax1.scatter(x2,y2)
for i in range(len(x1)):
ax1.text(x1[i],y1[i],"%s"%(i))
if i<len(x2):
ax1.text(x2[i],y2[i],"%s"%(i))

ax1.set_xlim(-3,3)
ax1.set_ylim(-3,3)
ax1.set_aspect("equal")

plt.show()

ana_aki.png

緑色の部分が反時計回りになっている所が重要。反時計回りはdx(差分)がマイナスになるので、負の面積というか負の断面二次モーメントが出てくる。

円形
N=21
x = np.cos(np.linspace(0, 2.0*np.pi, N))
y = np.sin(np.linspace(0, 2.0*np.pi, N))
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

maru.png

円は、理論値からの誤差がまだある。どのくらいサンプル間隔を細かくすればいいのか考えるために、下で誤差などを確認した。

円の理論値

$$\frac{d^{4}\pi}{64} \\ 半径r=1なので直径d=2\\ \frac{16\pi}{64}=0.25\pi \\ =‭0.78539816339744830961566084581988 \\ (電卓でπ=3.1415926535897932384626433832795‬を使い計算)‬$$

円周率にはnp.pi 3.141592653589793(16桁)を使うのでそれが原因の誤差が出るが、確認しない。

va = 0.25 * np.pi
Ix = []
Iy = []
Acx = []
Acy = []
Erx = []
Ery = []

for N in np.array(np.logspace(1,4,10),"int"):
#for N in np.array(np.logspace(1,4,10),"int")*4+1:
x = np.cos(np.linspace(0, 2.0*np.pi, N))
y = -np.sin(np.linspace(0, 2.0*np.pi, N))#Clockwise

I = simpson(x,y)
Ix.append(I[0])
Iy.append(I[1])

Acx.append(Ix[-1]/va)
Acy.append(Iy[-1]/va)

Erx.append(1 - Acx[-1])
Ery.append(1 - Acy[-1])
print("%s-points,\tIx:%s,\tIy:%s,\tAcx:%s,\tAcy:%s,\tErx:%s,\tEry:%s,\t"%(N,Ix[-1],Iy[-1],Acx[-1],Acy[-1],Erx[-1],Ery[-1]))
import matplotlib.pyplot as plt

x = np.array(np.logspace(1,4,10),"int")#*4+1
fig,(ax1,ax2) = plt.subplots(2,1)
ax1.plot(x,[va for i in range(len(x))])
ax1.scatter(x, Acy)

ax2.plot(x,[1 for i in range(len(x))])
ax2.scatter(x, Ery)

ax1.set_ylim(0.99,1)
ax2.set_ylim(0.00000000000001, 0.1)

ax1.set_xscale("log")
ax2.set_xscale("log")
ax1.set_yscale("log")
ax2.set_yscale("log")

plt.show()

figure_1 (3).png

上が精度で1に近づけばいい値。理論値で割って計算した。y軸はなんか知らんがすっ飛んだ。0.99~1.0間。

下が誤差で、0に近づけばいい値。1-精度で計算。両方とも対数グラフ。

サンプル数によっては、うまく計算できないのがある。実行して、printの結果を確認すればわかるけど、10の時に断面二次モーメントに異常な値が出る。dxの±が入れ替わる所か、yの±が入れ替わる所で起こっているのではないかと思っているが、確認は次回以降にする。とにかく、サンプル数に気を遣わずにシンプソン則を実行できなければならない。あと、欲を言えば、excelで出来れば言うことない。有料ソフトなのに、どこに行っても必ずあるし。このままでは、ベジェ曲線で自由に描いた線の積分ができない。円形で起こった現象は、ひし形でも出てきそう。正方形で出てこなかったのは、dxがぴったし0だから。

実用上は、40個位の分割でも円の断面二次モーメントを10-4の誤差で計算できる事が知れてよかった。

ちなみに、データ点数を4n+1(nは整数)にすれば、異常な値は出てこない。

現在の課題

  • ひし形、丸形などで、誤差が大きくなる。
  • 任意形状(ひらがなの「あ」とか)でやってない。

次回以降の目標

- 【断面二次】 シンプソン則のサンプル数を気にせず実行。

  • 【断面二次】Excelでやる。

  • 一覧表を作成し、誤差等をチェック
  • 【three.js】 pdbを表示する

絶対やるというわけではない。

posted by yuchan at 08:00 | Comment(51) | python

20190930

間に合わない。

断面二次モーメントと、数値積分のシンプソン則とかについて勉強したけど間に合わなかった。
three.jsについても勉強しているが、PDBLoaderというのが参考書通り動いてくれない。
月一では必ず書くつもりなので、まったく内容がなくても書く。
車の運転を練習し始めてからというものの、圧倒的に時間が足りなくなった。
いつもダラダラ勉強していると思っていたが、そんなことはなかったのか。。。

以上。
posted by yuchan at 07:00 | Comment(0) | 日記

20190902

three.jsのテスト

これは一記事にカウントしない。
断面二次(2)は、また別に、月末までにはやろうとしている。

スリージェイエスのテストをする。
あえてカタカナを割り振ったのは、世界に先駆けていると思う。
今後、3Dでぬるぬる動かしていく予定。three.jsのテスト

ローカルでテストしていたら、いい加減な書き方でも動いてしまった。
ブログにコードを貼り付けたら、そのまま動くと思いきや、動かなかった。

ブログのサービスも、javascriptをソースコードに書き込めないような仕様になっているのかもしれない。
コードを貼り付けてOKを押して、再度確認すると、変な文字が挿入されていた。
こういうのに気づくのにものすごい時間がかかる。変な文字入れるくらいなら、ポップアップ表示して、
「<script>タグは禁止されています」「<canvas>タグは挿入できません」とか書いてほしい。

どうしてもブログ内で動かしたい人は、codepenというサービスを使っている。
とにかく、このページでは動かない。
だから、外部リンクを貼った。


こういう黄色い棒がくるくる回るのを作ってみた。
ダウンロード (2).png
続きを読む
posted by yuchan at 19:00 | Comment(0) | three.js

20190831

断面二次モーメントを、座標点の配列から計算する(1)計算手順編

やること

断面二次モーメントを求める

断面二次モーメントを計算する方法は二つある。

断面形状の…

  1. 長さ寸法と理論式を使って計算する方法
  2. 座標値を使って、数値計算的に計算する方法。

がある。ここでは2をやる。

1と2の違いは、

  • 積分を形状別に解いて消してから、計算する場合

    と、

  • 積分を、狽ナ模擬してそのまま計算する場合

の違い。

ここでは、後者をやる。

任意形状の断面二次モーメントを求める

このページでは、閉じた線で描かれた任意形状の、x、y座標を使って、(面積や)断面二次モーメントを求める手順を考えてみた。

  1. 断面形状の縁を、閉じた線で描く様な、x,y座標を得る。
  2. 1を使って面積を求める。
  3. zに関する、もしくはyに関する断面一次モーメントを求める。
  4. 図心を求める。
    図心(y,z)=(z軸に関する断面一次モーメント, y軸に関する断面一次モーメント)/面積
  5. 1の座標から4の図心を引いた値を作り、その図形の断面二次モーメントを計算する。

普通は、断面二次モーメントの計算には理論式を使う。

断面二次モーメントは、普通、

  • 長さ寸法と理論式を使って計算する。
    これを実行するには…
    • 式を知っている必要がある。
      (→式が書いてある本等が要る)
      (下記amazon参照)
    • 断面形状によって式が違う。
    • 任意形状では式がない。
  • 理論式がない形状の場合、(solidworksなどの)有名なソフトの断面性能表示機能を使う。
  • 有償ソフトが使用できない場合、エクセルなどで自作する。(このページではexcelでやらない。pythonでやる。excelだって有料。)

のどれかで計算すると思う。
2つめ、3つめは、コンピュータで計算するので、理論値からの誤差がある。

←電子化(pdf(メニュー検索機能付き))全部入りDVD←

→1冊/全29冊 紙媒体 目的の式がある。→


↑断面二次モーメントの式が入っている。

長さ寸法を使うなら上の本などを使って計算されたし。


前述したが、閉じた線で描かれた任意形状の、x、y座標を使って、(面積や)断面二次モーメントを求める

1~5をやってみる

1. 文字データのx,y座標の取得

使用する座標値一覧は下の表。閉じた線を表す座標で、図形を表現している。これらの図形は、図心を中心にして作っていない。

図形 座標値
四角x[0,1,2,3,4,5,6,7,8,9,10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0]
 y[10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0,1,2,3,4,5,6,7,8,9,10]
せん断変換した四角x[41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0, 50.0, 51.0, 49.0, 47.0, 45.0, 43.0, 41.0, 39.0, 37.0, 35.0, 33.0, 31.0, 30.0, 29.0, 28.0, 27.0, 26.0, 25.0, 24.0, 23.0, 22.0, 21.0, 23.0, 25.0, 27.0, 29.0, 31.0, 33.0, 35.0, 37.0, 39.0, 41.0]
 y[17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 16.0, 15.0, 14.0, 13.0, 12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0]
45度回転した四角x[16.97056274847714, 17.67766952966369, 18.384776310850235, 19.09188309203678, 19.798989873223327, 20.506096654409877, 21.213203435596427, 21.920310216782973, 22.62741699796952, 23.334523779156065, 24.041630560342615, 23.334523779156065, 22.62741699796952, 21.920310216782973, 21.213203435596427, 20.506096654409877, 19.798989873223327, 19.09188309203678, 18.384776310850235, 17.67766952966369, 16.97056274847714, 16.26345596729059, 15.556349186104043, 14.849242404917495, 14.14213562373095, 13.435028842544401, 12.727922061357853, 12.020815280171306, 11.31370849898476, 10.606601717798211, 9.899494936611664, 10.606601717798211, 11.31370849898476, 12.020815280171306, 12.727922061357853, 13.435028842544401, 14.14213562373095, 14.849242404917495, 15.556349186104043, 16.26345596729059, 16.97056274847714]
 y[7.0710678118654755, 6.363961030678928, 5.65685424949238, 4.949747468305833, 4.242640687119286, 3.5355339059327378, 2.82842712474619, 2.1213203435596437, 1.4142135623730958, 0.7071067811865479, 0.0, -0.7071067811865479, -1.4142135623730958, -2.1213203435596437, -2.82842712474619, -3.5355339059327378, -4.242640687119286, -4.949747468305833, -5.65685424949238, -6.363961030678928, -7.0710678118654755, -6.363961030678928, -5.65685424949238, -4.949747468305832, -4.242640687119286, -3.5355339059327378, -2.82842712474619, -2.121320343559643, -1.4142135623730958, -0.7071067811865479, 0.0, 0.7071067811865479, 1.4142135623730958, 2.121320343559643, 2.82842712474619, 3.5355339059327378, 4.242640687119286, 4.949747468305832, 5.65685424949238, 6.363961030678928, 7.0710678118654755]
x[1.0, 0.94581724, 0.78914051, 0.54694816, 0.24548549, -0.08257935, -0.40169542, -0.67728157, -0.87947375, -0.9863613, -0.9863613, -0.87947375, -0.67728157, -0.40169542, -0.08257935, 0.24548549, 0.54694816, 0.78914051, 0.94581724, 1.0]
 y[0.0, 3.24699469e-01, 6.14212713e-01, 8.37166478e-01, 9.69400266e-01, 9.96584493e-01, 9.15773327e-01, 7.35723911e-01, 4.75947393e-01, 1.64594590e-01, -1.64594590e-01, -4.75947393e-01, -7.35723911e-01, -9.15773327e-01, -9.96584493e-01, -9.69400266e-01, -8.37166478e-01, -6.14212713e-01, -3.24699469e-01, 0.0]
穴あき四角x[[0,1,2,3,4,5,6,7,8,9,10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0],[1,1,1,1,1,1,1,1,1,2,3,4,5,6,7,8,9,9,9,9,9,9,9,9,9,8,7,6,5,4,3,2,1]]
 y[[10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0,1,2,3,4,5,6,7,8,9,10],[9,8,7,6,5,4,3,2,1,1,1,1,1,1,1,1,1,2,3,4,5,6,7,8,9,9,9,9,9,9,9,9,9]]
ひらがなの「あ」x[[1190, 1121.0, 1153, 1319.0, 1581.0, 1825.0, 1567.0, 1061.0, 1045, 1692.0, 1547.0, 1323.0, 1180.0, 910.0, 920.0, 930.0, 942.0, 899.0, 860.0, 797.0, 410.0, 244.0, 510.0, 797.0, 806.0, 811.0, 815, 717.0, 404.0, 443, 672.0, 824.0, 828, 830, 838.0, 738.0, 766, 981.0, 961.0, 930.0, 1182.0, 1186.0, 1243.0, 1321.0, 1405.0, 1434.0, 1393.0, 916.0, 893.0, 1190.0], [1194, 889.0, 887.0, 895.0, 1137.0, 1194.0], [793, 789.0, 790.0, 793.0, 621.0, 371.0, 459.0, 793.0]]
 y[[1008, 1094.0, 1130, 1006.0, 926.0, 541.0, 86.0, -74.0, -27, 541.0, 846.0, 938.0, 639.0, 332.0, 286.0, 250.0, 174.0, 57.0, 35.0, 240.0, 63.0, 268.0, 749.0, 913.0, 1054.0, 1124.0, 1176, 1171.0, 1333.0, 1366, 1251.0, 1262.0, 1307, 1343, 1462.0, 1618.0, 1661, 1522.0, 1421.0, 1280.0, 1356.0, 1358.0, 1393.0, 1419.0, 1399.0, 1366.0, 1321.0, 1186.0, 948.0, 1008.0], [942, 874.0, 680.0, 432.0, 786.0, 942.0], [334, 561.0, 663.0, 831.0, 719.0, 272.0, 163.0, 334.0]]
数字の「8」x[[908, 613.0, 740.0, 1043.0, 1322.0, 1459.0, 1143.0, 1528.0, 1401.0, 1028.0, 695.0, 535.0, 908.0], [1065, 1082.0, 1088, 1330.0, 1266.0, 1043.0, 877.0, 752.0, 850.0, 977.0, 1065.0], [977, 957.0, 668.0, 768.0, 1026.0, 1260.0, 1366.0, 1255.0, 1088.0, 977.0]]
 y[[856, 1225.0, 1477.0, 1573.0, 1495.0, 1231.0, 893.0, 457.0, 187.0, 74.0, 162.0, 453.0, 856.0], [924, 930.0, 934, 1231.0, 1413.0, 1507.0, 1460.0, 1229.0, 1043.0, 963.0, 924.0], [824, 814.0, 455.0, 234.0, 140.0, 222.0, 449.0, 664.0, 777.0, 824.0]]

上の表の座標値の作り方に関しては3種類あり…

  • 四角は、試行錯誤しながら手打ち。穴あきも苦労したが手打ち。四角だったから何とか打てた。蕎麦みたいな書きっぷり。
  • せん断、45度回転は、アフィン変換して作った。
  • 丸は、x座標はcosθ、y座標はsinθ (0<=θ<=2πを20分割(N=20))で作った。cosとかsinは、pythonのmathモジュールを使った。
  • ひらがなの「あ」、数字の「8」は、先月のページのやり方で作った。
    このページ下部に、フォントの文字データのx,y座標が取得のためのコードを載せました。
    (先月のと同じ内容です。IPAフォントのインストールを忘れずに。)

    フォントファイルのグリフのベクターを座標データに変換しています。

#!/usr/bin/env python
# -*- coding: cp932 -*-

x=ここに上表の右列のxを入力
y=ここに上表の右列のyを入力

こんな感じで変数x,yに入れておき、後のコードで使う。

変数x、yは、曲線のx座標y座標が入った配列です。
また、変数x2、y2の中には、一文字を表すために必要なパーツの分だけ、配列が入っています。
つまり、漢数字の「一」は一筆書できるので1つしか配列が入っていませんが、
同じ漢数字の「二」には配列が二つ(長さ2の配列)入っています。

上の座標を示して描ける線の特徴として、

  • 閉じた線
  • 外側が時計回り、内側が反時計回り

などがあります。

閉じた線

線を指定する座標は、最初の一点目の座標と同じ座標で終わるように作りました。線の途中で交差したりしていません。線の途中で交差すると、面積の計算に不都合があります。

時計回り、反時計回り。

閉じた線のx、y座標をあらわす配列が、時計回りに動いているか、それとも反時計回りに動いているかは、少なくとも下の二つの用途で重要です。

  • 文字の描画
    「の」や「○(まる)」などの穴開き文字(正式名称はあるのかすら知らない)に関しては、穴が開いている部分を表す配列が反時計回りになっている。下の図でいう緑の線。

の.png〇.png

  • 面積の計算
    ∫ydxすると、xが+に動くときdxが+に、xが-に動くときdxが-になる。これがうまい具合に効いて、閉じた曲線の面積が求まる。
    また、時計回り、反時計回りが切り替わると、負の面積が出ることがある。負の面積によって、文字の穴空き部分の面積が引き算される仕組みになっている。
    閉じた線の中で、交差があるとコレがうまくいかない。

断面形状の縁x、y座標をプロットする

2画以上の文字や、穴空き文字を表現するため、変数x,yは長さが1の配列ではない場合があります。

だから、matplotlib.pyplot.plotでそのまま使えません。matplotlib.pyplot.plotは、x,yの引数は、N*1のリストとかだからです。(Excelで折れ線グラフ一本書くときは、x,yそれぞれ一列ずつ選ぶのと一緒)

一つの変数の中に、たくさん作図すべき閉じた線が含まれる事があるので、

繰り返し構文+indexingで、plotしました。
上の表の座標を

x=

y=

とか書いた後に、下記コードの実行でx、yをプロットできます。

  • 必要なモジュール
    • matplotlib
    • datetime
    • time
#!/usr/bin/env python
# -*- coding: cp932 -*-

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'IPAGothic'

import time
from datetime import datetime as dt

f,(ax1) = plt.subplots(1,1)

if len(np.shape(y))==1:
ax1.plot(x,y)
ax1.scatter(x,y,facecolor=[0,0,0,0.3],edgecolor=[0,0,0])
else:
for i in range(len(x)):
ax1.plot(x[i],y[i])
ax1.scatter(x[i],y[i],facecolor=[0,0,0,0.3],edgecolor=[0,0,0])

ax1.set_aspect("equal")
#plt.show()

tdatetime = dt.now()
timtext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
fname = "%s.png"%(timtext)
plt.savefig(fname)

plt.close()

2.1を使って面積を求める。

x、yでプロットされた図形の面積は台形則で求める。

\int\int f(x,y) \ dA

面積を計算するだけなので、関数とか要らないけど、
f(x,y)=1

\int\int 1\ dydx

これを、台形則で模擬します。

#!/usr/bin/env python
# -*- coding: cp932 -*-

import numpy as np

A=0if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [(y[i]+y[i+1])/2.0 for i in range(len(y)-1)]
A = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
else:
for j in range(len(x)):#1文字に、複数の閉じた線が存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [(y[j][i]+y[j][i+1])/2.0 for i in range(len(y[j])-1)]
A1 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
A = A+A1

print A

曲線とは一致しないとは思いますが、曲線を間引いたデータとしては、正しい値が出てくると思います。もちろん細かく座標点をとればとるほど、正確な値に近づきます。
PCで扱うデータは、しょせん間が抜けた折れ線データに過ぎないので、面積に関しては 台形則で表してもいいと思います。

 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形square.pngsquare_sheared.pngsquare_45deg_rotated.png circle.pngsquare_hollow.pnggryph_a.pnggryph_8.png
面積100.0100.0100.0-3.08464495357…36694613.0402222.0

ついでに せん断変形で、面積が変わっていない事を確認

せん断写像のwikipediaにそう書いてあったきがする ので、せん断前後で面積に変化がない事を確認する。
せん断変換前後の図形のx,y座標をせん断変換(アフィン変換)で用意したので、
座標から面積を計算した結果を比較すると、せん断変形しようが、面積は変わらないことが確認できました。


もう理論値からズレてる。

曲線をPCで表現するとき、当然細かい折れ線で描く事になります。
間引いた折れ線データを取得した時点で、データに誤差が生じています。

だから、この後の操作がどんなに正しくとも、座標点として記録した時点で理論値からのズレはあります。
「それでもやるとしたら、どうすべきなの?」という観点で見て下さい。

半径1の円で、線の間引き方によって、値がずれていることが確認できるので載せました。

 粗い円やや粗い円細かい円理論値(半径2×3.14…)
図形circle_coarth2.pngcircle.pngcircle_precise.pngcircle_ideal.png
面積-2.89254424359…-3.08464495357…-3.14157194138…π≒3.14159265358…

3. zに関する、もしくはyに関する断面一次モーメントを求める。

y軸に関する断面一次モーメント $$\int y\ dA$$

x軸に関する断面一次モーメント$$\int x\ dA$$

断面一次モーメントは、図形がどの位置に描いてあるかによって、値が変わります。
原点から離れるほど大きくなります。

断面二次モーメントを計算する時も、図形が原点から離れるほど大きくなります。断面二次モーメントは、図心を原点と一致させてから求めることになります。次の手順4で、図形の重心(図心)を求めます。

断面一次モーメントは、その図心を求めるために計算します。

断面一次モーメントを求めるコードは下のように書いてみました。

#!/usr/bin/env python
# -*- coding: cp932 -*-
import numpy as np

C1=0
if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
C1 = C1+C
else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+1])/2.0) for i in range(len(y[j])-1)]
C2 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
C1 = C1+C2

print("X:%s"%(C1))

C3=0

if type(x[0])!=list:
dy = [y[i+1]-y[i] for i in range(len(x)-1)]
trapzx = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y)-1)])*-1
C3 = C3+C4

else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dy = [y[j][i+1]-y[j][i] for i in range(len(x[j])-1)]
trapzx = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+1])/2.0) for i in range(len(y[j])-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y[j])-1)])*-1
C3 = C3+C4

print("Y:%s"%(C3))
 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形square.pngsquare_sheared.pngsquare_45deg_rotated.pngcircle.pngsquare_hollow.pnggryph_a.pnggryph_8.png
断面一次モーメントX:500.0
Y:500.0
X:3600.0
Y:1200.0
X:1697.05627485
Y:1.42108547152e-14≒0
X:1.90819582357e-17≒0
Y:1.11022302463e-16≒0
X:180.0
Y:180.0
X:713082273.0
Y:475455483.5
X:419443127.5
Y:326005988.0

4. 図心を求める。

図心(x,y)=(y軸に関する断面一次モーメント, x軸に関する断面一次モーメント)/面積

あえて書く必要はない気がするけど、

  1. 先ず、面積を計算する。

  2. 断面一次モーメントを計算する

  3. 断面一次/面積

    してるのが下記。

#!/usr/bin/env python
# -*- coding: cp932 -*-

import
numpy as np
A=0

if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [(y[i]+y[i+1])/2.0 for i in range(len(y)-1)]
A = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
else:
for j in range(len(x)):#1文字に、複数の閉じた線が存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [(y[j][i]+y[j][i+1])/2.0 for i in range(len(y[j])-1)]
A1 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
A = A+A1
print A

C1=0

if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
C1 = C1+C
else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+ 1])/2.0) for i in range(len(y[j])-1)]
C2 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
C1 = C1+C2

print("X:%s"%(C1))

C3=0
if type(x[0])!=list:
dy = [y[i+1]-y[i] for i in range(len(x)-1)]
trapzx = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y)-1)])*-1
C3 = C3+C4
else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dy = [y[j][i+1]-y[j][i] for i in range(len(x[j])-1)]
trapzx = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+1])/2.0) for i in range(len(y[j])-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y[j])-1)])*-1
C3 = C3+C4

print("Y:%s"%(C3))
print
"First x moment of area%s"%(C1/A)
print "First y moment of area%s"%(C3/A)
 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形cen_square.pngcen_sheared_square.pngcen_45deg_rotated_square.pngcen_circle.pngcen_hollow_square.pngcen_ah.pngcen_8.png
断面一次モーメントX:500.0
Y:500.0
X:3600.0
Y:1200.0
X:1697.05627485
Y:1.42108547152e-14≒0
X:1.90819582357e-17
Y:1.11022302463e-16
X:180.0
Y:180.0
X:713082273.0
Y:475455483.5
X:419443127.5
Y:326005988.0
面積100.0100.0100.0-3.08464495357…36694613.0402222.0
図心X:5.0
Y:5.0
X:36.0
Y:12.0
X:16.9705627485
Y:1.42108547152e-16
X:6.18611169162e-18≒0
Y:3.59919225694e-17≒0
X:5.0
Y:5.0
X:1026.58929937
Y:684.489756886
X:1042.81498153
Y:810.512572659

5. 1の座標から4の図心を引いた値を作り、その図形の断面二次モーメントを計算する。

材料力学の教科書とか、ネット検索とから、
$$\int y^2\ dA$$

途中で $$\frac{1}{3} y^3$$ と、1/3が出てきたり、y3とか3乗になったりする点も気を付けたい。

積分の解き方に不安があったので、最初書いてたけど消した。

これをx,y座標から数値積分で求めたい。

先ず台形則で、積分を模擬してみた。(本当は二次の積分を台形則で模擬してはいけない気がしている。)


これまで、コードを載せてきたが、この断面二次モーメントの算出の仕方は、間違っているため、まだ載せない。

理論値からずれているのだ。

 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形cen_square.pngcen_sheared_square.pngcen_45deg_rotated_square.pngcen_circle.pngcen_hollow_square.pngcen_ah.pngcen_8.png
Ix=833.333333333
Iy=833.333333333
Ix=833.333333333
Iy=4233.33333333
Ix=841.666666667
Iy=841.666666667
Ix=0.771161239361
Iy=0.771161239361
Ix=492.0
Iy=492.0
Ix=1.21943806285e+11
Iy=1.026346451e+11
Ix=79177695468.9
Iy=38517544355.2
理論値からのズレなし-ありありなし--

(ずれていると判定するのに使った理論値は来月載せます。)

理論値からのずれ(誤差)は、なぜ起こるのか。積分を模擬する手法が台形則だからだと思う。

誤差がないものは、何でないのか?


水平、垂直に線が伸びている形状に関しては誤差がない。

斜めに線が伸びている形状は、誤差が大きい。


ここで、キレイなグラフを用意して説明しようとしたのだが、時間と通信残量制限との兼ね合いで、中途半端なままアップロードする。

これに関して改善するには、台形則(一次)で表現している部分をシンプソン則(二次)で表現する必要がありそうですが、来月に続きます。


今月使用した、断面二次モーメントの算出時の積分方法(台形則)と、

シンプソン則を比較したいと思います。多分、違ってくるはず、だと思いたい。不安になってきた。


今月は**「断面二次モーメントを、座標点の配列から計算する(1)計算手順編」**だったのに対し、

来月は**「断面二次モーメントを、座標点の配列から計算する(1)積分模擬変」**にしたいと思う。変って気づいたけど直さない。

8月中に間に合わんかったし、まだとちゅうだが、これは14000文字超えの超大作なので、まあゆるせ。というか間に合ってはいた。auの通信則で制限のせいだから。おらは悪くねぇだ。長続きしていると思う。

続きを読む
posted by yuchan at 23:59 | Comment(296) | python