fringeSearch.F90

Go to the documentation of this file.

subroutine fringesearch(oBlock, oFringe)
 
    use constants
    use block, only: fringetype
    use oversetdata, only: oversetblock, oversetfringe, localwallfringes, &
                           tmpfringeptr, nlocalwallfringe
    use inputoverset, onlY: overlapfactor, oversetprojtol
    use adtlocalsearch, only: mindistancetreesearchsinglepoint, &
                              containmenttreesearchsinglepoint
    use adtdata, only: adtbboxtargettype
    use adtutils, only: stack
    use utils, only: mynorm2
    use oversetutilities, only: fractoweights2, addtofringelist, addtofringebuffer, &
                                tic, toc, windindex
    implicit none
 
    type(oversetblock), intent(inout) :: oBlock
    type(oversetfringe), intent(inout) :: oFringe
 
    integer(kind=intType) :: idom, jdom
 
    ! Working Varaibles
    integer(kind=intType) :: nInterpol, elemID, nalloc, intInfo(3), intInfo2(3)
    integer(kind=intType) :: i, ii, jj, kk, j, nn, myI, myJ, myK
    integer(kind=intTYpe) :: iii, jjj, kkk, n, myind, nx, ny, nz, myindex
    logical :: invalid, failed
    real(kind=realtype) :: uu, vv, ww, err1, err2
    real(kind=realtype) :: uvw(5), uvw2(5), xx(4), pt(3), xcheck(3)
    real(kind=realtype), dimension(:, :), allocatable :: offset
    real(kind=realtype) :: oneminusu, oneminusv, oneminusw, weight(8)
    ! Variables we have to pass the ADT search routine
    integer(kind=intType), dimension(:), pointer :: BB
    type(adtbboxtargettype), dimension(:), pointer :: BB2
    integer(kind=intType), dimension(:), pointer :: frontLeaves
    integer(kind=intType), dimension(:), pointer :: frontLeavesNew
    type(fringetype) :: fringe
    ninterpol = 1 ! we get the ADT to compute the interpolated volume for us.
 
    ! Allocate the (pointer) memory that may be resized as necessary for
    ! the singlePoint search routine.
    allocate (bb(20), bb2(20), frontleaves(25), frontleavesnew(25), stack(100))
 
    ! Number of fringes we have:
    n = size(ofringe%x, 2)
 
    ! Offset vector:
    allocate (offset(3, n))
    offset = zero
    call tic(isurfacecorrection)
    call surfacecorrection(oblock, ofringe, offset, n)
    call toc(isurfacecorrection)
 
    call tic(idonorsearch)
    ! Search the cells one at a time:
    do i = 1, n
 
        ! Compute the potentailly offset point to search for.
        xx(1:3) = ofringe%x(:, i) + offset(:, i)
 
        call containmenttreesearchsinglepoint(oblock%ADT, xx, intinfo, uvw, &
                                              oblock%qualDonor, ninterpol, bb, frontleaves, frontleavesnew, failed)
 
        if (intinfo(1) >= 0) then
            call fractoweights2(uvw(1:3), weight)
            xcheck = zero
            do j = 1, 8
                xcheck = xcheck + weight(j) * oblock%xADT(:, oblock%hexaConn(j, intinfo(3)))
            end do
 
            if (mynorm2(xcheck - xx(1:3)) > oversetprojtol) then
                failed = .true.
            end if
        end if
 
        if (intinfo(1) >= 0 .and. failed) then
            ! we "found" a point but it is garbage. Do the failsafe search
            xx(4) = large
            call mindistancetreesearchsinglepoint(oblock%ADT, xx, intinfo, uvw, &
                                                  oblock%qualDonor, ninterpol, bb2, frontleaves, frontleavesnew)
 
            ! Check this one:
            call fractoweights2(uvw(1:3), weight)
            xcheck = zero
            do j = 1, 8
                xcheck = xcheck + weight(j) * oblock%xADT(:, oblock%hexaConn(j, intinfo(3)))
            end do
 
            ! Since this is the last line of defence, relax the tolerance a bit
            if (mynorm2(xcheck - xx(1:3)) > 100 * oversetprojtol) then
                ! This fringe has not found a donor
                intinfo(1) = -1
            else
                ! This one has now passed.
 
                ! Important! uvw(4) is the distance squared for this search
                ! not
                uvw(4) = uvw(5)
            end if
 
        end if
 
        elemfound: if (intinfo(1) >= 0) then
 
            ! Recompute the i,j,k indices on the donor
            elemid = intinfo(3) - 1 ! Make it zero based for the modding.
            ii = mod(elemid, oblock%il) + 1
            jj = mod(elemid / oblock%il, oblock%jl) + 1
            kk = elemid / (oblock%il * oblock%jl) + 1
 
            ! Now record the information onto the fringe
            fringe%donorProc = oblock%proc
            fringe%donorBlock = oblock%block
            fringe%dIndex = windindex(ii, jj, kk, oblock%il, oblock%jl, oblock%kl)
            fringe%donorFrac = uvw(1:3)
            fringe%quality = uvw(4)
 
            ! Also save the information about where it came from,
            ! we need this to combine everything together at the end.
            fringe%myBlock = ofringe%block
 
            myi = mod((i - 1), ofringe%nx) + 2
            myj = mod((i - 1) / ofringe%nx, ofringe%ny) + 2
            myk = (i - 1) / (ofringe%nx * ofringe%ny) + 2
            fringe%myIndex = windindex(myi, myj, myk, ofringe%il, ofringe%jl, ofringe%kl)
 
            ! Store the donor in the big flat list if it isn't invalid
            invalid = .false.
            do kkk = 0, 1
                do jjj = 0, 1
                    do iii = 0, 1
                        if (oblock%invalidDonor(ii + iii, jj + jjj, kk + kkk) .ne. 0) then
                            invalid = .true.
                        end if
                    end do
                end do
            end do
 
            if (.not. invalid) then
                call addtofringebuffer(ofringe%fringeIntBuffer, ofringe%fringeRealBuffer, &
                                       ofringe%nDonor, fringe)
            end if
 
            ! Save the fringe to the wallList. Note that we have to do
            ! this *after* the actual fringeList becuase we may modify the
            ! dI, dJ, dK here.
 
            if ((ofringe%isWall(i) > 0) .and. .not. (oblock%nearWall(ii, jj, kk) == 1)) then
 
                ! Here we have to recompute the i,j,k indices since we may
                ! need to modify them based on the frac.
 
                if (uvw(1) >= half) then
                    ii = ii + 1
                end if
 
                if (uvw(2) >= half) then
                    jj = jj + 1
                end if
 
                if (uvw(3) >= half) then
                    kk = kk + 1
                end if
 
                ! Recompute the full index for the wall.
                fringe%dIndex = windindex(ii, jj, kk, oblock%il, oblock%jl, oblock%kl)
 
                call addtofringelist(localwallfringes, nlocalwallfringe, fringe)
            end if
        end if elemfound
    end do
    deallocate (offset, bb, bb2, frontleaves, frontleavesnew, stack)
    call toc(idonorsearch)
end subroutine fringesearch