【发布时间】：2016-06-18 13:45:56
【问题描述】：
我有简单的代码，它用圆柱体包围的区域标记节点。在执行代码时，结果是圆柱体观察情况轻微倾斜 90 度
实际问题： 上述算法是在 Fortran 中实现的。如果在圆柱体内，代码会检查笛卡尔网格中的点。以下是测试用例： 圆柱在 yz 平面上相对于 y 轴成 90 度角。因此，方向向量$\vec{o}$为(0, 1, 0)。
案例 1： 方向向量直接用 $\vec{o}=(0.0,1.0,0.0)$ 赋值。这导致完美的圆柱体 $\theta=90.$
案例 2：
方向向量使用具有双精度精度的内在 Fortran 函数 dsin 和 dcos 指定 $\vec{o}=(0.0, \sin(\pi/2.0), \cos(\pi/2.0))$ $\pi$ 值分配有超过 20 个有效小数点。生成的圆柱体会导致轻微倾斜。
突出显示的区域表示由于圆柱体相对于笛卡尔轴倾斜而产生的额外材料。我也试过architecture specific maximum precision "pi" value.这也导致同样的问题。
这表明圆柱的实际角度不是90度。任何人都可以为这个问题提出有效的解决方案。我需要对任意角度使用内置三角函数并寻找准确的单元格标记方法。
注意：所有操作均以双精度精度执行。
实际功能如下。 rk 是定义参数，值为8
pure logical function in_particle(p,px,x)
type(md_particle_type),intent(in) :: p
real(kind=rk),intent(in) :: px(3),x(3)
real(kind=rk) :: r(3),rho(3),rop(2),ro2,rdiff,u
rop = particle_radii(p) ! (/R_orth,R_para/)
ro2 = rop(1)**2
rdiff = rop(2) - rop(1)
r = x-px
! Case 1:
! u = dot_product((/0.0_rk,-1.0_rk,0.0_rk/),r)
! rho = r-u*(/0.0_rk,-1.0_rk,0.0_rk/)
! Case 2:
u = dot_product((/0.0_rk,-dsin(pi/2.0_rk),dcos(pi/2.0_rk)/),r)
rho = r-u*(/0.0_rk,-dsin(pi/2.0_rk),dcos(pi/2.0_rk)/)
if((u.le.rdiff).and.(u.ge.-rdiff)) then
in_particle = dot_product(rho,rho) < ro2
else
in_particle = .false.
end if
end function in_particle
注意：三角运算是在代码内部完成的，以便更好地解释问题。然而，原始代码以矢量形式从用户那里读取方向。然后将此信息转换为四元数以进行粒子-粒子碰撞操作。在将四元数转换回方向向量时，这个错误甚至被放大了。甚至在碰撞开始之前，圆柱体的方向往往会被 2 个晶格单元所迷惑。
【问题讨论】：
-
请出示您的代码。
-
你的意思是说，如果你注释掉“Case 2”下面的两行，而取消注释“Case 1”下面的两行，你会得到预期的结果，而在上面的代码中你会得到一个意想不到的结果?如果插入
print *, dsin(pi/2.0_rk); stop之类的行（以验证 pi 的值是否正常）会得到什么？ -
只是一个风格点——dsin 和 dcos 非常适合 Fortran 66。简单的 sin 和 cos 完全没问题，已经有将近 40 年了。
-
再确认一下：你得到 cos(pi/2.0_rk) = 0.00000000000000006123... 吗？ （零的数量是 17。）如果是这样，我目前对这两行一无所知......至于其他行，我想知道 rdiff 是否保证为正（因为 IF 条件似乎假设它）。
-
虽然一个疯狂的猜测，给 IF 条件一些轻微的缓冲区，比如说，
- ( rdiff + eps ) < u .and. u < ( rdiff + eps )eps = 1.0d-3 等可能会有所帮助......？ （以避免可能的粒子恰好位于边界上的情况）。
标签： fortran geometry rounding computational-geometry