wallSearch.F90

Go to the documentation of this file.

module wallsearches
contains
    subroutine wallsearch(aSurf, bSurf)
 
        use constants
        use oversetdata, only: oversetwall, clusterareas
        use inputoverset, only: nearwalldist
        use adtlocalsearch, only: mindistancetreesearchsinglepoint
        use adtdata, only: adtbboxtargettype, adtleaftype
        use adtutils, only: stack
        use utils, only: mynorm2
        implicit none
 
        ! Input/Output
        type(oversetwall), intent(inout) :: aSurf, bSurf
 
        ! Working Varaibles
        integer(kind=intType) :: i, jj, k, iElem, maxLevels, nNeighbours, nOtherElem, iOther, otherElem
        integer(kind=intType) :: nInterpol, elemID, intInfo(3), factor, jelem, otherElems(4)
        real(kind=realtype) :: uvw(5), xx(4), dist, q1(3, 4), q2(3, 4), delta, radius1, radius2
 
        ! Variables we have to pass the ADT search routine
        integer(kind=intType), dimension(:), pointer :: frontLeaves
        integer(kind=intType), dimension(:), pointer :: frontLeavesNew
        type(adtbboxtargettype), dimension(:), pointer :: BB
        real(kind=realtype), dimension(3, 2) :: dummy
        integer(kind=intType), dimension(:), allocatable :: tmpCellArr
        integer(kind=intType), dimension(:), allocatable :: tmpNodeElem
        integer(kind=intType), dimension(:), allocatable :: mask
        type(adtleaftype), dimension(:), pointer :: ADTree
 
        logical :: overlapped
 
        if (asurf%nCells == 0 .or. bsurf%nCells == 0) then
            ! Either block doesn't have walls, so there is nothing do but just
            ! return.
            return
        end if
 
        if (clusterareas(bsurf%cluster) < clusterareas(asurf%cluster)) then
            ! B is smaller so we don't need to do anything
            return
        else if (clusterareas(bsurf%cluster) == clusterareas(asurf%cluster)) then
            ! If the areas are *exactly* equal (which can happen when we have
            ! replicated grids) cut by the cluster index.
 
            if (bsurf%cluster < asurf%cluster) then
                return
            end if
        end if
        ninterpol = 0
        ! Allocate the (pointer) memory that may be resized as necessary for
        ! the singlePoint search routine.
        allocate (bb(10), frontleaves(25), frontleavesnew(25), stack(100))
 
        ! Basically what we are doing it looping all of our bSurf NODES. We
        ! use a special "surface containment search". Essentially all we are
        ! looking for is if a point it inside of of an actual element
        ! BBox. If it isn't inside any BBox then we know it it can't
        ! overlap. This is essentialy fast cull of the majority of panels so
        ! we can later just focus on the ones that may actually overlap.
        ! node as being blanked
 
        ! Start with a max 10 layers (each with an unreduced 8 cells)
        maxlevels = 1
        allocate (tmpcellarr(3 * 3), mask(asurf%nCells), tmpnodeelem(bsurf%nNodes))
        tmpnodeelem(:) = 0
        mask = 0
        adtree => asurf%ADT%ADTree
 
        do i = 1, bsurf%nNodes
 
            xx(1:3) = bsurf%x(:, i)
            xx(4) = large
 
            ! Just check if it is inside the root bounding box..ie the full
            ! bounding box of the surface. This is would appear conservative,
            ! but isn't good enough. We need to expand by nearWallDist since
            ! it is possible a overlap occurs right at the edge of the
            ! bounding box.
            if (xx(1) >= adtree(1)%xMin(1) - nearwalldist .and. &
                xx(1) <= adtree(1)%xMax(4) + nearwalldist .and. &
                xx(2) >= adtree(1)%xMin(2) - nearwalldist .and. &
                xx(2) <= adtree(1)%xMax(5) + nearwalldist .and. &
                xx(3) >= adtree(1)%xMin(3) - nearwalldist .and. &
                xx(3) <= adtree(1)%xMax(6) + nearwalldist) then
 
                ! Now find the closest element on the other mesh for this
                ! node. This is the regular (expensive) closest point search
 
                call mindistancetreesearchsinglepoint(asurf%ADT, xx, intinfo, uvw, &
                                                      dummy, ninterpol, bb, frontleaves, frontleavesnew)
                elemid = intinfo(3)
                tmpnodeelem(i) = elemid
            end if
        end do
 
        ! Loop over the cells now since this is eventually want we need to blank out:
        cellloop: do i = 1, bsurf%nCells
 
            ! Extract out the elems found on the other mesh for the 4 nodes
            ! on my element. There could be none, 1 or up to 4 other elements.
 
            notherelem = 0
            do jj = 1, 4
                otherelem = tmpnodeelem(bsurf%conn(jj, i))
                if (otherelem /= 0) then
                    notherelem = notherelem + 1
                    otherelems(notherelem) = otherelem
                end if
            end do
 
            ! Get the coordinates of my quad
            do jj = 1, 4
                q1(:, jj) = bsurf%x(:, bsurf%conn(jj, i))
            end do
 
            do iother = 1, notherelem
 
                elemid = otherelems(iother)
 
                ! Get coordinates of the other (found) quad
                do jj = 1, 4
                    q2(:, jj) = asurf%x(:, asurf%conn(jj, elemid))
                end do
 
                ! Do a quick check of the cell itself. If it overlaps,
                ! we're done and don't need to deal with neighbor cells at
                ! all.
                call quadoverlap(q1, q2, overlapped)
                if (overlapped) then
                    bsurf%iBlank(bsurf%cellPtr(i)) = -2
 
                    ! No need to do anything else
                    cycle cellloop
                end if
 
                ! Otherwise, we need to do more work.
                radius1 = getcellradius(q1)
                radius2 = getcellradius(q2)
 
                ! We technically only should only need to add 1 here, but
                ! to be safer, we'll have at least two layers to check.
                factor = int(radius1 / radius2) + 2
 
                if (factor > maxlevels) then
                    deallocate (tmpcellarr)
                    maxlevels = factor
                    allocate (tmpcellarr((1 + 2 * maxlevels)**2))
                end if
 
                ! This is where it gets interesing: We can determine the
                ! number of recursive radiating layers we need to check
                ! based on the relative size of the the two quads.
 
                ! Now for the fun part: Recursion!
                nneighbours = 0
                call getneighbourcells(asurf, mask, elemid, factor, tmpcellarr, nneighbours)
 
                ! Now just blindly check them until we run out of find an overlapped one:
 
                do ielem = 1, nneighbours
                    elemid = tmpcellarr(ielem)
 
                    ! Return the mask for this elem back to 0
                    mask(elemid) = 0
 
                    ! Get coordinates of the other quad
                    do jj = 1, 4
                        q2(:, jj) = asurf%x(:, asurf%conn(jj, elemid))
                    end do
 
                    ! Do the actual overlap calc for the found cell:
                    call quadoverlap(q1, q2, overlapped)
 
                    if (overlapped) then
                        bsurf%iBlank(bsurf%cellPtr(i)) = -2
 
                        ! No need to do anything else, but we we do need to
                        ! flip all the mask elements back for the next
                        ! iteration of cellLoop
 
                        do jelem = ielem + 1, nneighbours
                            mask(tmpcellarr(jelem)) = 0
                        end do
 
                        cycle cellloop
 
                    end if
                end do
            end do
        end do cellloop
 
        deallocate (bb, frontleaves, frontleavesnew, stack, tmpcellarr, mask)
 
    end subroutine wallsearch
 
    recursive subroutine getneighbourcells(aSurf, mask, baseElemID, layers, elemList, nElemFound)
 
        ! This routine recursively assembles a list all neighbours within
        ! "layers" of the the baseELemID. The elemList is sorted such that
        ! there are no duplicates:
 
        use constants
        use oversetdata, only: oversetwall
        implicit none
 
        ! Input/Output
        type(oversetwall), intent(inout) :: asurf
        integer(kind=intType), intent(inout), dimension(:) :: mask, elemlist
        integer(kind=intType), intent(in) :: baseelemid, layers
        integer(kind=intType), intent(inout) :: nelemfound
 
        ! Working
        integer(kind=intType) :: i, inode, icell, curelem
 
        ! The recusive chain ends when layers == 0
        if (layers == 0) then
            return
        end if
 
        ! Loop over the nodes of the given quad:
        do i = 1, 4
            inode = asurf%conn(i, baseelemid)
 
            ! Loop over the (up to 4) cells surrounding this node use the
            ! node->elem (nte) array
            do icell = 1, 4
                curelem = asurf%nte(icell, inode)
                if (curelem /= 0) then
                    ! This is a real cell:
 
                    if (mask(curelem) /= baseelemid .and. mask(curelem) == 0) then
                        ! we know we don't need to add the baseElemID and if its
                        ! already in the mask we don't have to do anything either
 
                        nelemfound = nelemfound + 1
                        elemlist(nelemfound) = curelem
                        mask(curelem) = 1
                        ! Now recursively call again, with the baseElement of curElem and 1 fewer levels
                        call getneighbourcells(asurf, mask, curelem, layers - 1, elemlist, nelemfound)
                    end if
                end if
            end do
        end do
    end subroutine getneighbourcells
 
    subroutine quadoverlap(q1, q2, overlapped)
        ! Given two quad in *3D* determine if they overlap using the
        ! separation axis theorem after projecting onto the plane defined by
        ! the cell normal. Check both normals from each quad.
 
        use constants
        use utils, only: mynorm2, cross_prod
 
        implicit none
 
        ! input/output
        real(kind=realtype), dimension(3, 4), intent(in) :: q1, q2
        logical, intent(out) :: overlapped
 
        ! Working
        integer(kind=intType) :: ii, jj
        real(kind=realtype), dimension(2, 4) :: qq1, qq2
        real(kind=realtype), dimension(3) :: axis1, axis2, n1, n2, normal, v1, v2, c1, c2
        real(kind=realtype) :: e1, e2
        ! Check distance between cell centers
        c1 = zero
        c2 = zero
        do ii = 1, 4
            c1 = c1 + fourth * q1(:, ii)
            c2 = c2 + fourth * q2(:, ii)
        end do
 
        ! Get get max distance between center and node:
        e1 = zero
        e2 = zero
        do ii = 1, 4
            e1 = max(e1, mynorm2(c1 - q1(:, ii)))
            e2 = max(e2, mynorm2(c2 - q2(:, ii)))
        end do
 
        ! Check if distance between cell center sid beyond the threshold
        if (mynorm2(c1 - c2) .ge. (e1 + e2)) then
            overlapped = .false.
            return
        end if
 
        ! The two quads *may* be overlapped. We have to do it hard way.
 
        ! Normal of first quad
        v1 = q1(:, 3) - q1(:, 1)
        v2 = q1(:, 4) - q1(:, 2)
        call cross_prod(v1, v2, n1)
        n1 = n1 / mynorm2(n1)
 
        ! Normal of second quad
        v1 = q2(:, 3) - q2(:, 1)
        v2 = q2(:, 4) - q2(:, 2)
        call cross_prod(v1, v2, n2)
        n2 = n2 / mynorm2(n2)
 
        ! f the normals are not in the same direction, must be a thin
        ! surface.
        if (dot_product(n1, n2) < zero) then
            overlapped = .false.
            return
        end if
 
        do ii = 1, 2
            if (ii == 1) then
                normal = n1
                axis1 = q1(:, 2) - q1(:, 1)
            else
                normal = n2
                axis1 = q2(:, 2) - q2(:, 1)
            end if
 
            ! Project axis1 onto the plane and normalize
            axis1 = axis1 - dot_product(axis1, normal) * normal
            axis1 = axis1 / mynorm2(axis1)
 
            ! Axis 2 is now the normal cross axis1
            call cross_prod(normal, axis1, axis2)
            axis2 = axis2 / mynorm2(axis2)
 
            do jj = 1, 4
                qq1(1, jj) = dot_product(axis1, q1(:, jj))
                qq1(2, jj) = dot_product(axis2, q1(:, jj))
 
                qq2(1, jj) = dot_product(axis1, q2(:, jj))
                qq2(2, jj) = dot_product(axis2, q2(:, jj))
            end do
            call quadoverlap2d(qq1, qq2, overlapped)
 
            if (overlapped) then
                return
            end if
        end do
 
    end subroutine quadoverlap
 
    subroutine quadoverlap2d(q1, q2, overlapped)
        ! Given two quad in *2D* determine if they overlap using the
        ! separation axis theorem
 
        use constants
        implicit none
 
        ! input/output
        real(kind=realtype), dimension(2, 4), intent(in) :: q1, q2
        logical, intent(out) :: overlapped
 
        ! Working
        real(kind=realtype), dimension(4) :: tmp1, tmp2
        integer(kind=intType) :: ii, jj, kk, jjp1
        real(kind=realtype), dimension(2) :: axis, p0
        real(kind=realtype) :: min1, max1, min2, max2
        overlapped = .true.
        tmp1 = zero
        tmp2 = zero
        quadloop: do ii = 1, 2 ! Loop over the two quads
            edgeloop: do jj = 1, 4 ! Loop over the edges of each quad
                jjp1 = mod(jj, 4) + 1
 
                if (ii == 1) then
                    axis = q1(:, jjp1) - q1(:, jj)
                    p0 = q1(:, jj)
                else
                    axis = q2(:, jjp1) - q2(:, jj)
                    p0 = q2(:, jj)
                end if
 
                ! Take the axis normal
                axis = (/axis(2), -axis(1)/)
 
                ! Take the dot products
                do kk = 1, 4
                    tmp1(kk) = dot_product(axis, q1(:, kk) - p0)
                    tmp2(kk) = dot_product(axis, q2(:, kk) - p0)
                end do
 
                min1 = minval(tmp1)
                max1 = maxval(tmp1)
 
                min2 = minval(tmp2)
                max2 = maxval(tmp2)
 
                if (max1 < min2 .or. max2 < min1) then
                    overlapped = .false.
                    ! We can just jump right out since we know they cannot
                    ! overlap.
                    exit quadloop
                end if
            end do edgeloop
        end do quadloop
    end subroutine quadoverlap2d
 
    function getcellradius(q)
 
        use constants
        use utils, only: mynorm2
        implicit none
        ! Input
        real(kind=realtype), dimension(3, 4) :: q
        real(kind=realtype) :: getcellradius
 
        ! Working
        real(kind=realtype) :: c(3)
        integer(kind=intType) :: ii
 
        c = zero
        do ii = 1, 4
            c = c + fourth * q(:, ii)
        end do
 
        getcellradius = zero
        do ii = 1, 4
            getcellradius = max(getcellradius, mynorm2(c - q(:, ii)))
        end do
 
    end function getcellradius
 
end module wallsearches