Q&A
確かにその通りでした。gradientをとればいいのですものね。
気づかず恥ずかしいです。
それによって修正した場合，思ったとおりの画像を得られました。
最終的には任意の形状に対して法線ベクトルを計算できるようにしたいのですが，そのような場合，関数形がわからないのでgradientをとるわけにはいかないと思います。
そのような場合の解決方法をご存じでありませんか。
返答に質問で返して申し訳ありませんがよろしくお願いします
前提・実現したいこと
半球を生成し，その法線ベクトルを計算．
その後，与えられたV()というベクトルとその法線との積をとって，4Dmapとして表示したいと考えています．
発生している問題・エラーメッセージ
とりあえず計算は回っているようなのですが，得られる結果がおかしいです．
Vx=0,Vy=0,Vz=-100*7 なので，北極点で一番大きな値が得られると考えていたのですが，低い値が得られています．
法線ベクトルがうまく計算できていないのでしょうか．
x,y,z座標系での各格子点の距離の差をとって，ベクトルの成分とし，外積を計算して求めようとしているのですがこれでは上手くいかないのでしょうか．
また，半球が途中でぷつっと途切れているのも気になります．これは格子インデックスの最初と最後が一致していないからだと思うのですが，こちらも解答いただけると幸いです．
該当のソースコード
fortran
1 2 3program main 4 ! Determine some arrays 5 implicit none 6 integer,parameter :: i_max=100, j_max=100 7 double precision :: x(0:i_max,0:j_max), y(0:i_max,0:j_max), z(0:i_max,0:j_max) ! this is the position 8 double precision :: v(3) ! this is the velocity (uniform flow) 9 double precision :: cp(0:i_max,0:j_max) 10 integer :: i,j 11 12 call grid 13 call normal_and_Cp 14 call dataSave 15 print *, pi_generator() 16 stop 17 contains 18 19 double precision function pi_generator() 20 pi_generator = 4*atan(1.0) 21 return 22 end function 23 24 subroutine grid 25 ! Generating grids 26 double precision :: r = 2 ! radius 27 double precision :: theta(0:i_max) 28 double precision :: phi(0:j_max) 29 do 1000 j=0,j_max 30 theta(j)=float(j)/float((j_max))*(pi_generator() )/2.0 31 1000 continue 32 do 1001 i=0,i_max 33 phi(i) = float(i)/float(i_max)*(pi_generator() )*2 34 1001 continue 35 ! 球座標系で，r方向の変化はないモデルなので，2変数で書けそうだ 36 do 100 i=0,i_max 37 do 101 j=0,j_max 38 x(i,j) = r*cos(phi(i))*sin(theta(j)) 39 y(i,j) = r*sin(phi(i))*sin(theta(j)) 40 z(i,j) = r*cos(theta(j)) 41 101 continue 42 ! x(i,0)=x(i,j_max) 43 ! y(i,0)=y(i,j_max) 44 100 continue 45 end subroutine 46 47 subroutine normal_and_Cp 48 ! calculate cross product 49 ! 周り4格子点のx,y,zから法線ベクトルn(i,j)を求める． 50 ! まず，各格子点の周り4格子の座標の差をとり，ベクトルを作成する． 51 ! そのベクトルについて外積を計算．4つ分を足し合わせて，その平均をその格子が持つ法線ベクトルとする 52 ! 面と格子がずれることになるが，そこは一度妥協する 53 54 ! 55 ! o(U) 56 ! | 57 ! (L) | 58 ! o-------o-------o(R) 59 ! | 60 ! | 61 ! o(B) 62 63 double precision :: vec_R(3),vec_L(3), vec_U(3), vec_B(3), & 64 & cross_RU(3), cross_UL(3), cross_LB(3), cross_BR(3), & 65 & cross_avg(3) 66 do 200 i=0,i_max-1 67 do 201 j=0,j_max-1 68 ! i,j 格子点周りの成分を求めていく 69 ! 格子の作り方として，x,y座標を決めてからz座標をスイープ． 70 ! んで，x,y座標がthetaにそって更新されて，ほんでz座標スイープ 71 ! というようになっていることに注意 72 73 ! phi方向の変化を表すベクトルの成分 74 vec_R(1) = x(i+1,j) - x(i,j) ! x成分 75 vec_R(2) = y(i+1,j) - y(i,j) ! y成分 76 vec_R(3) = z(i+1,j) - z(i,j) ! z成分 77 ! print *, vec_R(3) ! phi方向なんで0.0でok 78 79 vec_L(1) = x(i-1,j) - x(i,j) 80 vec_L(2) = y(i-1,j) - y(i,j) 81 vec_L(3) = z(i-1,j) - z(i,j) 82 ! theta方向の変化を表すベクトルの成分 83 vec_U(1) = x(i,j+1) - x(i,j) 84 vec_U(2) = y(i,j+1) - y(i,j) 85 vec_U(3) = z(i,j+1) - z(i,j) 86 vec_B(1) = x(i,j-1) - x(i,j) 87 vec_B(2) = y(i,j-1) - y(i,j) 88 vec_B(3) = z(i,j-1) - z(i,j) 89 90 91 ! print *, vec_B(3) 92 ! 法線ベクトルは，R_Uの関係だと，成分は 93 cross_RU(1) = vec_R(2)*vec_U(3) - vec_R(3)*vec_U(2) 94 cross_RU(2) = vec_R(3)*vec_U(1)- vec_R(1)*vec_U(3) 95 cross_RU(3) = vec_R(1)*vec_U(2) - vec_R(2)*vec_U(1) 96 97 ! U_Lの関係だと，成分は， 98 cross_UL(1) = vec_U(2)*vec_L(3) - vec_U(3)*vec_L(2) 99 cross_UL(2) = vec_U(3)*vec_L(1) - vec_U(1)*vec_L(3) 100 cross_UL(3) = vec_U(1)*vec_L(2) - vec_U(2)*vec_L(1) 101 102 ! L_Bの関係だと，成分は， 103 cross_LB(1) = vec_L(2)*vec_B(3) - vec_L(3)*vec_B(2) 104 cross_LB(2) = vec_L(3)*vec_B(1) - vec_L(1)*vec_B(3) 105 cross_LB(3) = vec_L(1)*vec_B(2) - vec_L(2)*vec_B(1) 106 107 ! B_Rの関係だと，成分は， 108 cross_BR(1) = vec_B(2)*vec_R(3) - vec_B(3)*vec_R(2) 109 cross_BR(2) = vec_B(3)*vec_R(1) - vec_B(1)*vec_R(3) 110 cross_BR(3) = vec_B(1)*vec_R(2) - vec_B(2)*vec_R(1) 111 112 cross_avg(1) = (cross_RU(1) + cross_UL(1) + cross_LB(1) + cross_BR(1))/4.0 113 cross_avg(2) = (cross_RU(2) + cross_UL(2) + cross_LB(2) + cross_BR(2))/4.0 114 cross_avg(3) = (cross_RU(3) + cross_UL(3) + cross_LB(3) + cross_BR(3))/4.0 115 cross_avg(1) = cross_avg(1)/sqrt(cross_avg(1)*cross_avg(1)+1.0d-5) 116 cross_avg(2) = cross_avg(2)/sqrt(cross_avg(2)*cross_avg(2)+1.0d-5) 117 cross_avg(3) = cross_avg(3)/sqrt(cross_avg(3)*cross_avg(3)+1.0d-5) 118 119 ! 球の場合，normal=(x,y,z)でした．x**2+y**2+z**2=r**2より，f=x**2+y**2+z**2-r**2としてgradとればokでした... 120 ! ===================================================================================== 121 cross_avg(1)=0.5*x(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5) 122 cross_avg(2)=0.5*y(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5) 123 cross_avg(3)=0.5*z(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5) 124 ! ==================================================================================== 125 v(1)=0 126 v(2)=0 127 v(3)=-100*7 128 if (dot_product(v,cross_avg) < 0) then 129 cp(i,j) = 2*((dot_product(v,cross_avg)) & 130 & /sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)))**2 131 else 132 cp(i,j)=0 133 endif 134 201 continue 135 200 continue 136 end subroutine 137 138 double precision function dot_product(v,cross_avg) 139 double precision :: v(3), cross_avg(3) 140 dot_product=v(1)*cross_avg(1)+v(2)*cross_avg(2)+v(3)*cross_avg(3) 141 return 142 end function 143 144 subroutine dataSave 145 character (40):: filename 146 integer :: fo2=12 147 write(filename, "(a,i5.5,a)") "vtkdata/Newtonian" 148 open(fo2,file=filename) 149 do 700 i=1,i_max 150 do 701 j=1,j_max 151 write(fo2,*) x(i,j),y(i,j),z(i,j),cp(i,j) 152 701 continue 153 write(fo2,*) 154 700 continue 155 close(fo2) 156 end subroutine 157 158end program main
試したこと
北極点でメッシュがうまく存在しないのではと思い，格子数を増やしたり，角度をpi/2から変更したりしましたが上手くいきませんでした．
補足情報（FW/ツールのバージョンなど）
plotはgnuplotを使用して
splot 'Newtonian' u 1:2:3:4 w pm3d
で書いています．
fortran90を使用しています
2020/04/22 追記
法線を解析的に求めた場合の画像．
2020/04/23 追記
正規化方法をコメント頂いたように直した場合
回答2件
0
4つの外積を平均してから大きさで割っているように見えますが、4つの外積それぞれをそれぞれの大きさで割ってから平均しなければ駄目だと思います。
(訂正)単位化してから平均すると単位ベクトルでなくなるので、平均した後に更に単位化が必要でした。
また、どうやらそもそもこれが原因ではなかったようです。コメント参照。
投稿2020/04/22 15:24
編集2020/04/23 14:47
総合スコア3047
0
ベストアンサー
問題の現象と関係ないかもしれませんが……
法線ベクトルがうまく計算できていないのでしょうか．
x,y,z座標系での各格子点の距離の差をとって，ベクトルの成分とし，外積を計算して求めようとしているのですがこれでは上手くいかないのでしょうか．
何故そんな計算をせねばならないのかがわかりません．
球面上の位置(x,y,z)がそのまま法線なのではないのですか？
投稿2020/04/22 08:30
総合スコア11652
>任意の形状
がどう与えられるのか次第ではないでしょうか．
・形状が何らかの（微分可能な）連続関数として与えられるなら微分（導関数の式まではもらえない場合には数値微分になるのかな）から求められるかもしれません．
・そうではなくて，何か離散的なデータで与えられる場合（e.g. 点群データとか），近傍データから求めることになるのかな，と思います．（着目点の近傍に存在する座標点群に平面を当てはめる感じ？）
ありがとうございます．
点群データから求められるようにできればいいなと思っていました．
平面生成は最小二乗法等で平面を決めることになるでしょうか．
・まだ未定ですが，CADファイルを描きだして計算しようと思っています．
STLファイルで与えられた形状なら法線が既に与えられるようなので，うまくデータファイルをいじれるか挑戦してみようと思います．
個人の経験ですが，点群データに対して
> 最小二乗法
で法線を求めたことはあります．法線を計算したい着目位置から近傍各点までの距離に応じて適当な重みづけした最小二乗で．
データが綺麗じゃない（何かのセンサでの実測値だとかいう）場合には対ノイズ性を考慮した工夫は必要だとは思いますが，人工データならまぁ大抵はそんなんで行けるんじゃないかと．
CAD由来だとデータがそもそも（三角形とかの）「面の集合」になると思うので法線方向の算出で苦労することは無いかな，と．
（３ＤＣＧの世界だと「法線を補間」だとかいう別の話も有り得ますが）
STLファイルでやってみたところ，それらしい結果を得ることができそうです．gnuplotでSTLファイルをどうプロットするかというところにまだ難がありますが，法線データが求まってるものを使うのは確実性が高そうだなと感じています．
・教えていただいた最小二乗法，今の自分にはまだ難しそうですが．やってみたいと思います．
> 最小二乗法
単純に座標群に平面を当てはめるなら「結果平面のパラメータとは～なんやかんやで～最小固有値に対応する固有ベクトルである」とかいう話になり，「固有値固有ベクトルを求めてくれるライブラリか何かがあればまぁ実装いける」という感じになると思います．
（２次元の直線当てはめの例とかを検索すればこんな感じの話がたくさん見つかるでしょう．）
「俺はそんな話は一切無視して数値探索でもそもそとやってやるぜ！」という方向も（処理速度とか結果精度とかが許すならば）有りかもです．（←やってみるだけなら，割と楽しい）
あと，gnuplotが当てはめやってくれる機能を持っていた気もします．
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2020/04/23 01:26
2020/04/23 01:34
2020/04/23 01:38
2020/04/23 13:23
2020/04/23 14:39
2020/04/23 14:40
2020/04/24 01:16
2020/04/24 01:25 編集
2020/04/27 01:13
2020/04/27 01:32 編集
2020/04/30 13:24