前提・実現したいこと
半球を生成し,その法線ベクトルを計算.
その後,与えられた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を使用しています
回答2件
あなたの回答
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