|
IGMBaseLib 1.0
|
src/util/math/FVM_HDiff_operator.f90
Go to the documentation of this file.
00001 00011 module FVM_HDiff_operator 00012 00013 ! モジュール引用 ; Use statement 00014 ! 00015 00016 ! 種類型パラメータ 00017 ! Kind type parameter 00018 ! 00019 use dc_types, only: DP ! 倍精度実数型. Double precision. 00020 00021 ! 線形代数 00022 ! Linear algebra 00023 ! 00024 use igmcore_linear_algebra, only: & 00025 & vec_length, vec_normarize, & 00026 & cross, dot 00027 00028 ! 球面三角法 00029 ! Spherical trigonometry 00030 ! 00031 use igmcore_spherical_trigonometry, only: & 00032 & geodesic_arc_length 00033 00034 ! 宣言文 ; Declaration statement 00035 ! 00036 implicit none 00037 private 00038 00039 ! 公開手続き 00040 ! Public procedure 00041 ! 00042 public :: eval_div_operator, eval_rot_operator, eval_grad_operator 00043 00044 contains 00045 00046 00047 ! 00051 function eval_div_operator( & 00052 & val_GP0, val_CV, GP0, CV, & ! (in) 00053 & CVArea, ic_radius, CV_num, val_dim, ret_dim & ! (in) 00054 & ) result(val) 00055 00056 ! 宣言文 ; Declaration statements 00057 ! 00058 00059 integer, intent(in) :: val_dim 00060 integer, intent(in) :: ret_dim 00061 integer, intent(in) :: CV_num 00062 real(DP), intent(in) :: val_GP0(val_dim) 00063 real(DP), intent(in) :: GP0(3) 00064 real(DP), intent(in) :: val_CV(CV_num, val_dim) 00065 real(DP), intent(in) :: CV(CV_num, val_dim) 00066 real(DP), intent(in) :: CVArea 00067 real(DP), intent(in) :: ic_radius 00068 real(DP) val(ret_dim) 00069 00070 ! 実行文 ; Executable statements 00071 ! 00072 00073 val(:) = div_operator( & 00074 & vec_CV= val_CV, pt_CV=CV, CVArea=CVArea, & ! (in) 00075 & pt_size=CV_num, ic_radius=ic_radius ) ! (in) 00076 00077 end function eval_div_operator 00078 00079 ! 00083 function eval_rot_operator( & 00084 & val_GP0, val_CV, GP0, CV, & ! (in) 00085 & CVArea, ic_radius, CV_num, val_dim, ret_dim & ! (in) 00086 & ) result(val) 00087 00088 ! 宣言文 ; Declaration statements 00089 ! 00090 integer, intent(in) :: val_dim 00091 integer, intent(in) :: ret_dim 00092 integer, intent(in) :: CV_num 00093 real(DP), intent(in) :: val_GP0(val_dim) 00094 real(DP), intent(in) :: GP0(3) 00095 real(DP), intent(in) :: val_CV(CV_num, val_dim) 00096 real(DP), intent(in) :: CV(CV_num, val_dim) 00097 real(DP), intent(in) :: CVArea 00098 real(DP), intent(in) :: ic_radius 00099 real(DP) val(ret_dim) 00100 00101 ! 実行文 ; Executable statements 00102 ! 00103 00104 val(:) = rot_operator( & 00105 & vec_CV= val_CV, pt_CV=CV, CVArea=CVArea, & ! (in) 00106 & pt_size=CV_num, ic_radius=ic_radius ) ! (in) 00107 00108 end function eval_rot_operator 00109 00110 ! 00114 function eval_grad_operator( & 00115 & val_GP0, val_CV, GP0, CV, & ! (in) 00116 & CVArea, ic_radius, CV_num, val_dim, ret_dim & ! (in) 00117 & ) result(val) 00118 00119 ! 宣言文 ; Declaration statements 00120 ! 00121 integer, intent(in) :: val_dim 00122 integer, intent(in) :: ret_dim 00123 integer, intent(in) :: CV_num 00124 real(DP), intent(in) :: val_GP0(val_dim) 00125 real(DP), intent(in) :: GP0(3) 00126 real(DP), intent(in) :: val_CV(CV_num, val_dim) 00127 real(DP), intent(in) :: CV(CV_num, val_dim) 00128 real(DP), intent(in) :: CVArea 00129 real(DP), intent(in) :: ic_radius 00130 real(DP) val(ret_dim) 00131 00132 ! 実行文 ; Executable statements 00133 ! 00134 00135 ! val(:) = grad_operator( & 00136 ! & q_GP0=val_GP0(1), q_CV= val_CV(:,1), pt_CV=CV, CVArea=CVArea, pt_size=CV_num, ic_radius=ic_radius ) 00137 00138 00139 val(:) = grad_operator2( & 00140 & q_GP0=val_GP0(1), q_CV= val_CV(:,1), GP=GP0, pt_CV=CV, & ! (in) 00141 & CVArea=CVArea, pt_size=CV_num, ic_radius=ic_radius ) ! (in) 00142 00143 end function eval_grad_operator 00144 00145 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00146 ! 00147 ! 非公開手続き 00148 ! 00149 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00150 00151 ! 00161 function div_operator( vec_CV, pt_CV, CVArea, pt_size, ic_radius ) result(val) 00162 00163 ! 宣言文 ; Declaration statement 00164 ! 00165 integer, intent(in) :: pt_size 00166 real(DP), intent(in) :: vec_CV(pt_size, 3) 00167 real(DP), intent(in) :: pt_CV(pt_size, 3) 00168 real(DP), intent(in) :: CVArea 00169 real(DP), intent(in) :: ic_radius 00170 real(DP) val 00171 00172 ! 作業変数 00173 ! Work variables 00174 ! 00175 integer :: i 00176 real(DP) :: n(3) 00177 00178 ! 実行文 ; Executable statement 00179 ! 00180 00181 val = 0.0d0 00182 00183 do i=1, pt_size 00184 n(:) = cross( pt_CV(1+mod(i, pt_size), :), pt_CV(i,:) ) 00185 call vec_normarize(n) 00186 00187 val = val + geodesic_arc_length(ic_radius, pt_CV(1+mod(i, pt_size), :), pt_CV(i,:)) * & 00188 & dot( 0.5d0 * ( vec_CV(1+mod(i, pt_size), :) + vec_CV(i,:) ), n(:) ) 00189 end do 00190 00191 val = val / CVArea 00192 00193 end function div_operator 00194 00195 ! 00205 function rot_operator( vec_CV, pt_CV, CVArea, pt_size, ic_radius ) result(val) 00206 00207 ! 宣言文 ; Declaration statement 00208 ! 00209 integer, intent(in) :: pt_size 00210 real(DP), intent(in) :: vec_CV(pt_size, 3) 00211 real(DP), intent(in) :: pt_CV(pt_size, 3) 00212 real(DP), intent(in) :: CVArea 00213 real(DP), intent(in) :: ic_radius 00214 real(DP) val 00215 00216 ! 作業変数 00217 ! Work variable 00218 ! 00219 integer :: i 00220 real(DP) :: m(3) 00221 00222 ! 実行文 ; Executable statement 00223 ! 00224 val = 0.0d0 00225 00226 do i=1, pt_size 00227 m(:) = pt_CV(1+mod(i, pt_size), :) - pt_CV(i,:) 00228 call vec_normarize(m) 00229 00230 val = val + geodesic_arc_length(ic_radius, pt_CV(1+mod(i, pt_size), :), pt_CV(i,:)) * & 00231 & dot( 0.5d0 * ( vec_CV(1+mod(i, pt_size), :) + vec_CV(i,:) ), m(:) ) 00232 end do 00233 00234 val = val / CVArea 00235 00236 end function rot_operator 00237 00238 ! 00249 function grad_operator( q_GP0, q_CV, pt_CV, CVArea, pt_size, ic_radius ) result(val) 00250 00251 ! 宣言文 ; Declaration statement 00252 ! 00253 integer, intent(in) :: pt_size 00254 real(DP), intent(in) :: q_GP0 00255 real(DP), intent(in) :: q_CV(pt_size) 00256 real(DP), intent(in) :: pt_CV(pt_size, 3) 00257 real(DP), intent(in) :: CVArea 00258 real(DP), intent(in) :: ic_radius 00259 real(DP) val(3) 00260 00261 ! 作業変数 00262 ! Work variable 00263 ! 00264 integer :: i 00265 real(DP) :: n(3) 00266 00267 ! 実行文 ; Executable statement 00268 ! 00269 00270 val(:) = 0.0d0 00271 00272 do i=1, pt_size 00273 n(:) = cross( pt_CV(1+mod(i, pt_size), :), pt_CV(i,:) ) 00274 call vec_normarize(n) 00275 00276 val(:) = val(:) + geodesic_arc_length(ic_radius, pt_CV(1+mod(i, pt_size), :), pt_CV(i,:)) * & 00277 & ( 0.5d0 * ( q_CV(1+mod(i, pt_size)) + q_CV(i) ) - q_GP0 ) * n(:) 00278 end do 00279 00280 val(:) = val(:) / CVArea 00281 00282 end function grad_operator 00283 00284 ! 00295 function grad_operator2( q_GP0, q_CV, GP, pt_CV, CVArea, pt_size, ic_radius ) result(val) 00296 00297 ! 宣言文 ; Declaration statement 00298 ! 00299 integer, intent(in) :: pt_size 00300 real(DP), intent(in) :: q_GP0 00301 real(DP), intent(in) :: q_CV(pt_size) 00302 real(DP), intent(in) :: GP(3) 00303 real(DP), intent(in) :: pt_CV(pt_size, 3) 00304 real(DP), intent(in) :: CVArea 00305 real(DP), intent(in) :: ic_radius 00306 real(DP) val(3) 00307 00308 ! 作業変数 00309 ! Work variables 00310 ! 00311 integer :: i 00312 real(DP) :: n(3) 00313 00314 ! 実行文 ; Executable statement 00315 ! 00316 00317 val(:) = 0.0d0 00318 00319 do i=1, pt_size 00320 n(:) = cross( pt_CV(1+mod(i, pt_size), :) - pt_CV(i,:), GP(:) ) 00321 call vec_normarize(n) 00322 00323 val(:) = val(:) + geodesic_arc_length(ic_radius, pt_CV(1+mod(i, pt_size), :), pt_CV(i,:)) * & 00324 & ( 0.5d0 * ( q_CV(1+mod(i, pt_size)) + q_CV(i) ) - q_GP0 ) * n(:) 00325 end do 00326 00327 val(:) = val(:) / CVArea 00328 00329 end function grad_operator2 00330 00331 end module FVM_HDiff_operator
