いわて駐在研究日誌2。

NEVER STAND BEHIND ME

VF計算(その8)

また放置状態。スラドのほうは日記を消してしまったので、ここで書く。
VF計算で、合い向かうメッシュ同士の間に、他のメッシュがあるかどうかを判定するプログラムを考え中だが、
あんまり速くない....。どうしたものか。

    do j=1, n_bndfc
    do i=1, j-1
      if(isfaced(i,j,1) == 1) then

        ! make bounding box of i and j
        xmin = min(xyzb(i,1,1),xyzb(i,2,1),xyzb(i,3,1),xyzb(j,1,1),xyzb(j,2,1),xyzb(j,3,1))
        xmax = max(xyzb(i,1,1),xyzb(i,2,1),xyzb(i,3,1),xyzb(j,1,1),xyzb(j,2,1),xyzb(j,3,1))

        ymin = min(xyzb(i,1,2),xyzb(i,2,2),xyzb(i,3,2),xyzb(j,1,2),xyzb(j,2,2),xyzb(j,3,2))
        ymax = max(xyzb(i,1,2),xyzb(i,2,2),xyzb(i,3,2),xyzb(j,1,2),xyzb(j,2,2),xyzb(j,3,2))

        zmin = min(xyzb(i,1,3),xyzb(i,2,3),xyzb(i,3,3),xyzb(j,1,3),xyzb(j,2,3),xyzb(j,3,3))
        zmax = max(xyzb(i,1,3),xyzb(i,2,3),xyzb(i,3,3),xyzb(j,1,3),xyzb(j,2,3),xyzb(j,3,3))

        ! make ray from i to j
        ray (1:3) = cgb(j,1:3)-cgb(i,1:3)
        rayn(1:3) = normalize(ray(1:3))
        
        do k=1, n_bndfc
          if(k /= i .and. k /= j                     .and. &
             cgb(k,1) >= xmin .and. cgb(k,1) <= xmax .and. &
             cgb(k,2) >= ymin .and. cgb(k,2) <= ymax .and. &
             cgb(k,3) >= zmin .and. cgb(k,3) <= zmax       )then
            
            wkvec1(1:3) = cgb(k,1:3) - cgb(i,1:3)
            wkvec2(1:3) = -  wkvec1(1:3) + dot_product(rayn(1:3),wkvec1(1:3))*rayn(1:3)
            rk2 = dot_product(wkvec2(1:3),wkvec2(1:3))  ! square of normal distance of cgb(k) from ray(i->j)

            if( rk2 <= edgemax2(k) ) then              
              call is_raycross(cgb(i,1:3),rayn(1:3),xyzb(k,1:3,1:3),cgb(k,1:3),nvecb(k,1:3),edgebvec(k,1:3,1:3),isobs)
              if(isobs > 0) then
                isobstacle(i,j) = isobs
                isobstacle(j,i) = isobs
                exit
              endif
            endif

          endif
        enddo
         
      endif
    enddo
    enddo