adtBuild.F90

Go to the documentation of this file.

module adtbuild
    !
    !      Module which contains all the subroutines for the building of
    !      an ADT, both surface and volume.
    !
    use constants
    use adtutils, only: adtleaftype, adts, stack, nstack, adtterminate, &
                        qsortbboxes, reallocateadts, allocateadts
    use adtdata, only: adttype
    implicit none
 
    !=================================================================
 
contains
 
    !===============================================================
 
    subroutine buildadt(ADT)
        !
        !        This routine builds the 6 dimensional ADT for the given
        !        ADT. When this routine is called it is assumed that the
        !        bounding boxes of the grid have already been computed; the
        !        ADT for these bounding boxes is built here.
        !        Subroutine intent(inout) arguments.
        !        --------------------------------
        !        ADT: adt derived type to build
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        type(adttype), intent(inout) :: ADT
        !
        !       Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: i, j, k, ii, kk, ll, mm, nn, nfl, nfr
        integer(kind=intType) :: nLeaves, nBBoxes, splitDir
        integer(kind=intType) :: nLeavesToDivide, nLeavesToDivideNew
        integer(kind=intType) :: nLeavesTot
 
        integer(kind=intType), dimension(:), pointer :: BB_IDs
        integer(kind=intType), dimension(:), pointer :: BB_IDsNew
        integer(kind=intType), dimension(:), pointer :: nBB_IDs
        integer(kind=intType), dimension(:), pointer :: nBB_IDsNew
        integer(kind=intType), dimension(:), pointer :: curLeaf
        integer(kind=intType), dimension(:), pointer :: curLeafNew
        integer(kind=intType), dimension(:), pointer :: tmpIntPointer
 
        integer(kind=intType), dimension(0:ADT%nProcs - 1) :: tmpArr
 
        real(kind=realtype), dimension(:, :), pointer :: xbbox
 
        real(kind=realtype), dimension(3, 2) :: rootleafbbox
        real(kind=realtype), dimension(3, 2, 0:ADT%nProcs - 1) :: &
            rootleavesbbox
 
        type(adtleaftype), dimension(:), pointer :: adtree
 
        ! Initialize nStack and allocate the corresponding array stack.
        ! These are used in the qsort routine for the bounding boxes.
        ! As this routine is called quite often it is more efficient to
        ! have the stack array available rather than allocate it over
        ! and over again.
 
        nstack = 100
        allocate (stack(nstack), stat=ierr)
        if (ierr /= 0) &
            call adtterminate(adt, "buildADT", &
                              "Memory allocation failure for stack.")
 
        ! Determine the number of leaves of the adt. It can be proved
        ! that nLeaves equals nBBoxes - 1 for an optimally balanced
        ! tree. Take the exceptional case of nBBoxes == 0 and
        ! nBBoxes == 1 into account.
 
        nbboxes = adt%nBBoxes
        nleaves = nbboxes - 1
        if (nbboxes <= 1) nleaves = nleaves + 1
 
        adt%nLeaves = nleaves
 
        ! Allocate the memory for the adt.
 
        allocate (adt%ADTree(nleaves), stat=ierr)
        if (ierr /= 0) &
            call adtterminate(adt, "buildADT", &
                              "Memory allocation failure for ADTree.")
 
        ! Set some pointers to make the code more readable.
 
        xbbox => adt%xBBox
        adtree => adt%ADTree
 
        ! Allocate the memory for the arrays which control the
        ! subdivision of the leaves.
 
        nn = (nbboxes + 1) / 2
        nn = max(nn, 1_inttype)
 
        allocate (bb_ids(nbboxes), bb_idsnew(nbboxes), &
                  nbb_ids(0:nn), nbb_idsnew(0:nn), &
                  curleaf(nn), curleafnew(nn), stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "buildADT", &
             "Memory allocation failure for the arrays &
             &used in the subdivision.")
 
        ! Initialize the arrays BB_IDs, nBB_IDs and curLeaf, such that
        ! all bounding boxes belong to the root leaf. Also set the
        ! counters nLeavesToDivide and nLeavesTot, depending on the
        ! situation
 
        nbb_ids(0) = 0; nbb_ids(1) = nbboxes
        curleaf(1) = 1
 
        do i = 1, nbboxes
            bb_ids(i) = i
        end do
 
        nleavestodivide = min(nleaves, 1_inttype)
        nleavestot = nleavestodivide
 
        ! Initialize splitDir to 0, such that the first time it will
        ! split in direction 1.
 
        splitdir = 0
 
        ! Loop to subdivide the leaves. The division is such that the
        ! adt is optimally balanced.
 
        leafdivision: do
 
            ! Criterion to exit the loop.
 
            if (nleavestodivide == 0) exit
 
            ! Initializations for the next round of subdivisions and
            ! increment splitDir.
 
            nleavestodividenew = 0
            nbb_idsnew(0) = 0
 
            splitdir = splitdir + 1
            if (splitdir > 6) splitdir = 1
 
            ! Loop over the current number of leaves to be divided.
 
            currentleavesloop: do i = 1, nleavestodivide
 
                ! Store the number of bounding boxes present in the leaf
                ! in nn, the current leaf number in mm and i-1 in ii.
 
                ii = i - 1
                nn = nbb_ids(i) - nbb_ids(ii)
                mm = curleaf(i)
 
                ! Determine the bounding box coordinates of this leaf.
 
                ll = bb_ids(nbb_ids(ii) + 1)
                adtree(mm)%xMin(1) = xbbox(1, ll)
                adtree(mm)%xMin(2) = xbbox(2, ll)
                adtree(mm)%xMin(3) = xbbox(3, ll)
                adtree(mm)%xMin(4) = xbbox(4, ll)
                adtree(mm)%xMin(5) = xbbox(5, ll)
                adtree(mm)%xMin(6) = xbbox(6, ll)
 
                adtree(mm)%xMax(1) = xbbox(1, ll)
                adtree(mm)%xMax(2) = xbbox(2, ll)
                adtree(mm)%xMax(3) = xbbox(3, ll)
                adtree(mm)%xMax(4) = xbbox(4, ll)
                adtree(mm)%xMax(5) = xbbox(5, ll)
                adtree(mm)%xMax(6) = xbbox(6, ll)
 
                do j = (nbb_ids(ii) + 2), nbb_ids(i)
                    ll = bb_ids(j)
 
                    adtree(mm)%xMin(1) = min(adtree(mm)%xMin(1), xbbox(1, ll))
                    adtree(mm)%xMin(2) = min(adtree(mm)%xMin(2), xbbox(2, ll))
                    adtree(mm)%xMin(3) = min(adtree(mm)%xMin(3), xbbox(3, ll))
                    adtree(mm)%xMin(4) = min(adtree(mm)%xMin(4), xbbox(4, ll))
                    adtree(mm)%xMin(5) = min(adtree(mm)%xMin(5), xbbox(5, ll))
                    adtree(mm)%xMin(6) = min(adtree(mm)%xMin(6), xbbox(6, ll))
 
                    adtree(mm)%xMax(1) = max(adtree(mm)%xMax(1), xbbox(1, ll))
                    adtree(mm)%xMax(2) = max(adtree(mm)%xMax(2), xbbox(2, ll))
                    adtree(mm)%xMax(3) = max(adtree(mm)%xMax(3), xbbox(3, ll))
                    adtree(mm)%xMax(4) = max(adtree(mm)%xMax(4), xbbox(4, ll))
                    adtree(mm)%xMax(5) = max(adtree(mm)%xMax(5), xbbox(5, ll))
                    adtree(mm)%xMax(6) = max(adtree(mm)%xMax(6), xbbox(6, ll))
                end do
 
                ! Determine the situation of the leaf. It is either a
                ! terminal leaf or a leaf that must be subdivided.
 
                terminaltest: if (nn <= 2) then
 
                    ! Terminal leaf. Store the ID's of the bounding boxes with
                    ! negative numbers in children.
 
                    adtree(mm)%children(1) = -bb_ids(nbb_ids(ii) + 1)
                    adtree(mm)%children(2) = -bb_ids(nbb_ids(i))
 
                    else terminaltest
 
                    ! Leaf must be divided. Sort the bounding boxes of the
                    ! current leaf in increasing order; the sorting is based
                    ! on the coordinate in the split direction.
 
                    call qsortbboxes(bb_ids(nbb_ids(ii) + 1:), nn, adt, splitdir)
 
                    ! Determine the number of bounding boxes in the left leaf.
                    ! This number is at least 2. The actual number stored in
                    ! kk is this number plus an offset. Also initialize the
                    ! counter nfl, which is used to store the bounding boxes
                    ! in the arrays for the new round.
 
                    kk = (nn + 1) / 2 + nbb_ids(ii)
                    nfl = nbb_idsnew(nleavestodividenew)
 
                    ! Copy the ID's of the left bounding boxes into BB_IDsNew.
                    ! Also update nLeavesToDivideNew and the corresponding
                    ! entry in nBB_IDsNew.
 
                    do k = (nbb_ids(ii) + 1), kk
                        nfl = nfl + 1
                        bb_idsnew(nfl) = bb_ids(k)
                    end do
 
                    nleavestodividenew = nleavestodividenew + 1
                    nbb_idsnew(nleavestodividenew) = nfl
 
                    ! Update the total number of leaves and store this number
                    ! in child 1 of the current leaf and in the current leaves
                    ! for the next round.
 
                    nleavestot = nleavestot + 1
                    adtree(mm)%children(1) = nleavestot
                    curleafnew(nleavestodividenew) = nleavestot
 
                    ! The right leaf will only be created if it has more than
                    ! one bounding box in it, i.e. if the original leaf has
                    ! more than three bounding boxes. If the new leaf only has
                    ! one bounding box in it, it is not created; instead the
                    ! bounding box is stored in the current leaf.
 
                    if (nn == 3) then
 
                        ! Only three bounding boxes present in the current leaf.
                        ! The right leaf is not created and the last bounding
                        ! box is stored as the second child of the current leaf.
 
                        adtree(mm)%children(2) = -bb_ids(nbb_ids(i))
 
                    else
 
                        ! More than 3 bounding boxes are present and thus the
                        ! right leaf is created. Copy the ID's from BB_IDs into
                        ! BB_IDsNew and update the counters for the new round.
 
                        nfr = nbb_idsnew(nleavestodividenew)
                        do k = (kk + 1), nbb_ids(i)
                            nfr = nfr + 1
                            bb_idsnew(nfr) = bb_ids(k)
                        end do
 
                        nleavestodividenew = nleavestodividenew + 1
                        nbb_idsnew(nleavestodividenew) = nfr
 
                        ! Update the total number of leaves and store this number
                        ! in child 2 of the current leaf and in the current
                        ! leaves for the next round.
 
                        nleavestot = nleavestot + 1
                        adtree(mm)%children(2) = nleavestot
                        curleafnew(nleavestodividenew) = nleavestot
 
                    end if
 
                end if terminaltest
 
            end do currentleavesloop
 
            ! Swap the pointers for the next round.
 
            nleavestodivide = nleavestodividenew
 
            tmpintpointer => bb_ids
            bb_ids => bb_idsnew
            bb_idsnew => tmpintpointer
 
            tmpintpointer => nbb_ids
            nbb_ids => nbb_idsnew
            nbb_idsnew => tmpintpointer
 
            tmpintpointer => curleaf
            curleaf => curleafnew
            curleafnew => tmpintpointer
 
        end do leafdivision
 
        ! Deallocate the arrays used to build the local tree.
 
        deallocate (stack, bb_ids, bb_idsnew, nbb_ids, nbb_idsnew, &
                    curleaf, curleafnew, stat=ierr)
        if (ierr /= 0) &
            call adtterminate(adt, "buildADT", &
                              "Deallocation failure for the local arrays.")
        !
        !        Local tree has been built. Now determine the global
        !        information for this tree.
        !
        ! Determine the number and processor ID's of the non-empty local
        ! trees by gathering the data from the participating processors.
 
        ii = 1
        if (nbboxes == 0) ii = 0
 
        call mpi_allgather(ii, 1, adflow_integer, tmparr, 1, adflow_integer, &
                           adt%comm, ierr)
 
        ii = 0
        do i = 0, (adt%nProcs - 1)
            ii = ii + tmparr(i)
        end do
 
        adt%nRootLeaves = ii
 
        ! Allocate the memory for both the processor ID's and the
        ! 3D bounding box of the root leaf.
 
        allocate (adt%rootLeavesProcs(ii), &
                  adt%rootBBoxes(3, 2, ii), stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "buildADT", &
             "Memory allocation failure for &
             &rootLeavesProcs and rootBBoxes.")
 
        ! Determine the processor ID's of the non-empty trees.
 
        ii = 0
        adt%myEntryInRootProcs = 0
 
        do i = 0, (adt%nProcs - 1)
            if (tmparr(i) > 0) then
                ii = ii + 1
                adt%rootLeavesProcs(ii) = i
                if (adt%myID == i) adt%myEntryInRootProcs = ii
            end if
        end do
 
        ! Determine the local 3D bounding box of the root leaf.
        ! If no local tree is present, just set it to zero to avoid
        ! problems. The values do not matter.
 
        if (nbboxes == 0) then
            rootleafbbox = zero
        else
            rootleafbbox(1, 1) = adtree(1)%xMin(1)
            rootleafbbox(2, 1) = adtree(1)%xMin(2)
            rootleafbbox(3, 1) = adtree(1)%xMin(3)
 
            rootleafbbox(1, 2) = adtree(1)%xMax(4)
            rootleafbbox(2, 2) = adtree(1)%xMax(5)
            rootleafbbox(3, 2) = adtree(1)%xMax(6)
        end if
 
        ! Gather the data of the root leaves.
 
        call mpi_allgather(rootleafbbox, 6, adflow_real, rootleavesbbox, &
                           6, adflow_real, adt%comm, ierr)
 
        ! Store the 3D root bounding boxes of the non-empty trees in
        ! the data structure for the current ADT.
 
        ii = 0
        do i = 0, (adt%nProcs - 1)
            if (tmparr(i) > 0) then
                ii = ii + 1
 
                adt%rootBBoxes(1, 1, ii) = rootleavesbbox(1, 1, i)
                adt%rootBBoxes(2, 1, ii) = rootleavesbbox(2, 1, i)
                adt%rootBBoxes(3, 1, ii) = rootleavesbbox(3, 1, i)
 
                adt%rootBBoxes(1, 2, ii) = rootleavesbbox(1, 2, i)
                adt%rootBBoxes(2, 2, ii) = rootleavesbbox(2, 2, i)
                adt%rootBBoxes(3, 2, ii) = rootleavesbbox(3, 2, i)
            end if
        end do
 
    end subroutine buildadt
 
    subroutine buildsurfaceadt(nTria, nQuads, nNodes, &
                               coor, triaConn, quadsConn, &
                               BBox, useBBox, comm, &
                               adtID)
        !
        !        This routine builds the 6 dimensional ADT, which stores the
        !        given surface grid. The memory intensive part of these
        !        arguments, the arrays with the coordinates and
        !        connectivities, are not copied. Instead pointers are set to
        !        these arrays. It is therefore the responsibility of the user
        !        not to deallocate this memory before all the searches have
        !        been performed.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nNodes:    Number of local nodes in the given grid.
        !        nTria:     Idem for the triangles.
        !        nQuads:    Idem for the quadrilaterals.
        !        BBox(3,2): The possible bounding box. Only elements within
        !                   this box will be stored in the ADT.
        !        useBBox:   Whether or not to use the bounding box.
        !        comm:      MPI-communicator for the global ADT.
        !        adtID:     The ID of the ADT.
        !        Subroutine intent(in), target arguments.
        !        ----------------------------------------
        !        coor(3,nNodes):      Nodal coordinates of the local grid.
        !        triaConn(3,nTria):   Local connectivity of the triangles.
        !        quadsConn(4,nQuads): Idem for the quadrilaterals.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer, intent(in) :: comm
        character(len=*), intent(in) :: adtid
 
        integer(kind=intType), intent(in) :: ntria
        integer(kind=intType), intent(in) :: nquads
        integer(kind=intType), intent(in) :: nnodes
 
        logical, intent(in) :: usebbox
 
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: triaconn
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: quadsconn
 
        real(kind=realtype), dimension(3, 2), intent(in) :: bbox
 
        real(kind=realtype), dimension(:, :), intent(in), &
            target :: coor
        !
        !       Local variables.
        !
        integer :: ierr, ll, nnpe, jj
 
        integer(kind=adtElementType) :: eltype
 
        integer(kind=intType) :: i, j, ii, mm, nn, nelem
 
        integer(kind=intType), dimension(:, :), pointer :: conn
 
        real(kind=realtype), dimension(3) :: xmin, xmax
 
        logical, dimension(:), allocatable :: elementwithinbbox
 
        type(adttype), pointer :: ADT
 
        ! ! Allocate or reallocate the memory for ADTs. This depends
        ! ! whether or not this is the first ADT to be built.
 
        if (allocated(adts)) then
            call reallocateadts(adtid, jj)
        else
            call allocateadts
            jj = 1
        end if
 
        ! The ADT we are currently working with
        adt => adts(jj)
 
        ! Make sure the ADT is active and store the ID of this ADT.
 
        adt%isActive = .true.
        adt%adtID = adtid
 
        ! Copy the communicator and determine the number of processors
        ! and my processor ID in this group.
 
        adt%comm = comm
        call mpi_comm_rank(comm, adt%myID, ierr)
        call mpi_comm_size(comm, adt%nProcs, ierr)
 
        ! Set the ADT type, which is a surface ADT.
 
        adt%adtType = adtsurfaceadt
 
        ! Copy the number of nodes and surface elements and set the
        ! number of volume elements to 0; only a surface grid has been
        ! given.
 
        adt%nNodes = nnodes
        adt%nTria = ntria
        adt%nQuads = nquads
 
        adt%nTetra = 0
        adt%nPyra = 0
        adt%nPrisms = 0
        adt%nHexa = 0
 
        ! Set the pointers for the coordinates and the
        ! surface connectivities.
 
        adt%coor => coor
        adt%triaConn => triaconn
        adt%quadsConn => quadsconn
 
        ! Determine the number of elements to be stored in the ADT.
        ! This depends whether or not the global bounding box should be
        ! used when building the ADT.
 
        testbbox: if (usebbox) then
 
            ! Global bounding box is used. Allocate the memory for the
            ! logical elementWithinBBox.
 
            nn = ntria + nquads
            allocate (elementwithinbbox(nn), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildSurfaceADT", &
                 "Memory allocation failure for &
                 &elementWithinBBox.")
 
            ! Loop over the number of element types.
 
            ii = 0
            elementloop1: do ll = 1, 2
 
                ! Set the correct pointers for this element.
 
                call setsurfacepointers(ll)
 
                ! Loop over the elements and determine the bounding box of
                ! each element.
 
                do i = 1, nelem
                    ii = ii + 1
 
                    mm = conn(1, i)
                    xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
                    xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
                    xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
                    do j = 2, nnpe
                        mm = conn(j, i)
 
                        xmin(1) = min(xmin(1), coor(1, mm))
                        xmin(2) = min(xmin(2), coor(2, mm))
                        xmin(3) = min(xmin(3), coor(3, mm))
 
                        xmax(1) = max(xmax(1), coor(1, mm))
                        xmax(2) = max(xmax(2), coor(2, mm))
                        xmax(3) = max(xmax(3), coor(3, mm))
                    end do
 
                    ! Check if the bounding box is (partially) inside the
                    ! global bounding box. If so, set elementWithinBBox
                    ! to .true.; otherwise set it to .false.
 
                    if (xmax(1) >= bbox(1, 1) .and. xmin(1) <= bbox(1, 2) .and. &
                        xmax(2) >= bbox(2, 1) .and. xmin(2) <= bbox(2, 2) .and. &
                        xmax(3) >= bbox(3, 1) .and. xmin(3) <= bbox(3, 2)) then
                        elementwithinbbox(ii) = .true.
                    else
                        elementwithinbbox(ii) = .false.
                    end if
 
                end do
            end do elementloop1
 
            ! Determine the local number of elements within the global
            ! bounding box.
 
            ii = 0
            do i = 1, nn
                if (elementwithinbbox(i)) ii = ii + 1
            end do
 
            adt%nBBoxes = ii
 
            ! Allocate the memory for the bounding box coordinates, the
            ! corresponding element type and the index in the connectivity.
 
            allocate (adt%xBBox(6, ii), adt%elementType(ii), &
                      adt%elementID(ii), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildSurfaceADT", &
                 "Memory allocation failure for bounding &
                 &box data.")
 
            ! Repeat the loop over the all the elements, but now store
            ! the bounding boxes in the ADT.
 
            ii = 0
            nn = 0
            elementloop2: do ll = 1, 2
 
                ! Set the correct pointers for this element.
 
                call setsurfacepointers(ll)
 
                ! Loop over the elements and store the bounding box info,
                ! if needed.
 
                do i = 1, nelem
                    ii = ii + 1
                    testwithin: if (elementwithinbbox(ii)) then
 
                        nn = nn + 1
 
                        adt%elementType(nn) = eltype
                        adt%elementID(nn) = i
 
                        mm = conn(1, i)
                        xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
                        xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
                        xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
                        do j = 2, nnpe
                            mm = conn(j, i)
 
                            xmin(1) = min(xmin(1), coor(1, mm))
                            xmin(2) = min(xmin(2), coor(2, mm))
                            xmin(3) = min(xmin(3), coor(3, mm))
 
                            xmax(1) = max(xmax(1), coor(1, mm))
                            xmax(2) = max(xmax(2), coor(2, mm))
                            xmax(3) = max(xmax(3), coor(3, mm))
                        end do
 
                        adt%xBBox(1, nn) = xmin(1)
                        adt%xBBox(2, nn) = xmin(2)
                        adt%xBBox(3, nn) = xmin(3)
 
                        adt%xBBox(4, nn) = xmax(1)
                        adt%xBBox(5, nn) = xmax(2)
                        adt%xBBox(6, nn) = xmax(3)
 
                    end if testwithin
                end do
            end do elementloop2
 
            ! Deallocate the memory for elementWithinBBox.
 
            deallocate (elementwithinbbox, stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildSurfaceADT", &
                 "Deallocation failure for &
                 &elementWithinBBox.")
 
            else testbbox
 
            ! No global bounding box. The number of local bounding boxes
            ! to be stored is the total number of local surface elements.
 
            ii = ntria + nquads
            adt%nBBoxes = ii
 
            ! Allocate the memory for the bounding box coordinates, the
            ! corresponding element type and the index in the connectivity.
 
            allocate (adt%xBBox(6, ii), adt%elementType(ii), &
                      adt%elementID(ii), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildSurfaceADT", &
                 "Memory allocation failure for bounding &
                 &box data.")
 
            ! Loop over the number of element types present, i.e. 2,
            ! to store the bounding boxes; nn is the counter.
 
            nn = 0
            elementloop3: do ll = 1, 2
 
                ! Set the correct pointers for this element.
 
                call setsurfacepointers(ll)
 
                ! Loop over the number of elements and store the bounding
                ! box info.
 
                do i = 1, nelem
                    nn = nn + 1
 
                    adt%elementType(nn) = eltype
                    adt%elementID(nn) = i
 
                    mm = conn(1, i)
                    xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
                    xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
                    xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
                    do j = 2, nnpe
                        mm = conn(j, i)
 
                        xmin(1) = min(xmin(1), coor(1, mm))
                        xmin(2) = min(xmin(2), coor(2, mm))
                        xmin(3) = min(xmin(3), coor(3, mm))
 
                        xmax(1) = max(xmax(1), coor(1, mm))
                        xmax(2) = max(xmax(2), coor(2, mm))
                        xmax(3) = max(xmax(3), coor(3, mm))
                    end do
 
                    adt%xBBox(1, nn) = xmin(1)
                    adt%xBBox(2, nn) = xmin(2)
                    adt%xBBox(3, nn) = xmin(3)
 
                    adt%xBBox(4, nn) = xmax(1)
                    adt%xBBox(5, nn) = xmax(2)
                    adt%xBBox(6, nn) = xmax(3)
                end do
 
            end do elementloop3
 
        end if testbbox
 
        ! Build the ADT from the now known boundary boxes.
 
        call buildadt(adt)
 
        !===============================================================
 
    contains
 
        !=============================================================
 
        subroutine setsurfacepointers(ll)
            !
            !          This internal subroutine sets the pointers to the correct
            !          surface element, such that a loop over the element types
            !          can be used.
            !          Subroutine intent(in) arguments.
            !          --------------------------------
            !          ll: Element type for which the pointers must be used.
            !
            implicit none
            !
            !         Subroutine arguments.
            !
            integer, intent(in) :: ll
 
            select case (ll)
            case (1)
                eltype = adttriangle; nelem = ntria; nnpe = 3
                conn => triaconn
            case (2)
                eltype = adtquadrilateral; nelem = nquads; nnpe = 4
                conn => quadsconn
            case (3)
            end select
 
        end subroutine setsurfacepointers
 
    end subroutine buildsurfaceadt
 
    subroutine buildvolumeadt(nTetra, nPyra, nPrisms, &
                              nHexa, nNodes, coor, &
                              tetraConn, pyraConn, prismsConn, &
                              hexaConn, BBox, useBBox, &
                              comm, adtID)
        !
        !        This routine builds the 6 dimensional ADT, which stores the
        !        given volume grid. The memory intensive part of these
        !        arguments, the arrays with the coordinates and
        !        connectivities, are not copied. Instead pointers are set to
        !        these arrays. It is therefore the responsibility of the user
        !        not to deallocate this memory before all the searches have
        !        been performed.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nNodes:    Number of local nodes in the given grid.
        !        nTetra:    Idem for the tetrahedra.
        !        nPyra:     Idem for the pyramids.
        !        nPrisms:   Idem for the prisms.
        !        nHexa:     Idem for the hexahedra.
        !        BBox(3,2): The possible bounding box. Only elements within
        !                   this box will be stored in the ADT.
        !        useBBox:   Whether or not to use the bounding box.
        !        comm:      MPI-communicator for the global ADT.
        !        adtID:     The ID of the ADT.
        !        Subroutine intent(in), target arguments.
        !        ----------------------------------------
        !        coor(3,nNodes):        Nodal coordinates of the local grid.
        !        tetraConn(4,nTetra):   Local connectivity of the tetrahedra.
        !        pyraConn(5,nPyra):     Idem for the pyramids.
        !        prismsConn(6,nPrisms): Idem for the prisms.
        !        hexaConn(8,nHexa):     Idem for the hexahedra.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer, intent(in) :: comm
        character(len=*), intent(in) :: adtID
 
        integer(kind=intType), intent(in) :: nTetra
        integer(kind=intType), intent(in) :: nPyra
        integer(kind=intType), intent(in) :: nPrisms
        integer(kind=intType), intent(in) :: nHexa
        integer(kind=intType), intent(in) :: nNodes
 
        logical, intent(in) :: useBBox
 
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: tetraConn
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: pyraConn
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: prismsConn
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: hexaConn
 
        real(kind=realtype), dimension(3, 2), intent(in) :: bbox
 
        real(kind=realtype), dimension(:, :), intent(in), &
            target :: coor
        !
        !       Local variables.
        !
        integer :: ierr, ll, nNPE
 
        integer(kind=adtElementType) :: elType
 
        integer(kind=intType) :: i, j, ii, jj, mm, nn, nElem
 
        integer(kind=intType), dimension(:, :), pointer :: conn
 
        real(kind=realtype), dimension(3) :: xmin, xmax
 
        logical, dimension(:), allocatable :: elementWithinBBox
 
        type(adttype), pointer :: ADT
 
        ! Allocate or reallocate the memory for ADTs. This depends
        ! whether or not this is the first ADT to be built.
 
        if (allocated(adts)) then
            call reallocateadts(adtid, jj)
        else
            call allocateadts
            jj = 1
        end if
 
        adt => adts(jj)
 
        ! Make sure the ADT is active and store the ID of this ADT.
 
        adt%isActive = .true.
        adt%adtID = adtid
 
        ! Copy the communicator and determine the number of processors
        ! and my processor ID in this group.
 
        adt%comm = comm
        call mpi_comm_rank(comm, adt%myID, ierr)
        call mpi_comm_size(comm, adt%nProcs, ierr)
 
        ! Set the ADT type, which is a volume ADT.
 
        adt%adtType = adtvolumeadt
 
        ! Copy the number of nodes and volume elements and set the number
        ! of surface elements to 0; only a volume grid has been given.
 
        adt%nNodes = nnodes
        adt%nTetra = ntetra
        adt%nPyra = npyra
        adt%nPrisms = nprisms
        adt%nHexa = nhexa
 
        adt%nTria = 0
        adt%nQuads = 0
 
        ! Set the pointers for the coordinates and the
        ! volume connectivities.
 
        adt%coor => coor
        adt%tetraConn => tetraconn
        adt%pyraConn => pyraconn
        adt%prismsConn => prismsconn
        adt%hexaConn => hexaconn
 
        ! Determine the number of elements to be stored in the ADT.
        ! This depends whether or not the global bounding box should be
        ! used when building the ADT.
 
        testbbox: if (usebbox) then
 
            ! Global bounding box is used. Allocate the memory for the
            ! logical elementWithinBBox.
 
            nn = ntetra + npyra + nprisms + nhexa
            allocate (elementwithinbbox(nn), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildVolumeADT", &
                 "Memory allocation failure for &
                 &elementWithinBBox.")
 
            ! Loop over the number of element types.
 
            ii = 0
            elementloop1: do ll = 1, 4
 
                ! Set the correct pointers for this element.
 
                call setvolumepointers(ll)
 
                ! Loop over the elements and determine the bounding box of
                ! each element.
 
                do i = 1, nelem
                    ii = ii + 1
 
                    mm = conn(1, i)
                    xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
                    xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
                    xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
                    do j = 2, nnpe
                        mm = conn(j, i)
 
                        xmin(1) = min(xmin(1), coor(1, mm))
                        xmin(2) = min(xmin(2), coor(2, mm))
                        xmin(3) = min(xmin(3), coor(3, mm))
 
                        xmax(1) = max(xmax(1), coor(1, mm))
                        xmax(2) = max(xmax(2), coor(2, mm))
                        xmax(3) = max(xmax(3), coor(3, mm))
                    end do
 
                    ! Check if the bounding box is (partially) inside the
                    ! global bounding box. If so, set elementWithinBBox
                    ! to .true.; otherwise set it to .false.
 
                    if (xmax(1) >= bbox(1, 1) .and. xmin(1) <= bbox(1, 2) .and. &
                        xmax(2) >= bbox(2, 1) .and. xmin(2) <= bbox(2, 2) .and. &
                        xmax(3) >= bbox(3, 1) .and. xmin(3) <= bbox(3, 2)) then
                        elementwithinbbox(ii) = .true.
                    else
                        elementwithinbbox(ii) = .false.
                    end if
 
                end do
            end do elementloop1
 
            ! Determine the local number of elements within the global
            ! bounding box.
 
            ii = 0
            do i = 1, nn
                if (elementwithinbbox(i)) ii = ii + 1
            end do
 
            adt%nBBoxes = ii
 
            ! Allocate the memory for the bounding box coordinates, the
            ! corresponding element type and the index in the connectivity.
 
            allocate (adt%xBBox(6, ii), adt%elementType(ii), &
                      adt%elementID(ii), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildVolumeADT", &
                 "Memory allocation failure for bounding &
                 &box data.")
 
            ! Repeat the loop over the all the elements, but now store
            ! the bounding boxes in the ADT.
 
            ii = 0
            nn = 0
            elementloop2: do ll = 1, 4
 
                ! Set the correct pointers for this element.
 
                call setvolumepointers(ll)
 
                ! Loop over the elements and store the bounding box info,
                ! if needed.
 
                do i = 1, nelem
                    ii = ii + 1
                    testwithin: if (elementwithinbbox(ii)) then
 
                        nn = nn + 1
 
                        adt%elementType(nn) = eltype
                        adt%elementID(nn) = i
 
                        mm = conn(1, i)
                        xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
                        xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
                        xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
                        do j = 2, nnpe
                            mm = conn(j, i)
 
                            xmin(1) = min(xmin(1), coor(1, mm))
                            xmin(2) = min(xmin(2), coor(2, mm))
                            xmin(3) = min(xmin(3), coor(3, mm))
 
                            xmax(1) = max(xmax(1), coor(1, mm))
                            xmax(2) = max(xmax(2), coor(2, mm))
                            xmax(3) = max(xmax(3), coor(3, mm))
                        end do
 
                        adt%xBBox(1, nn) = xmin(1)
                        adt%xBBox(2, nn) = xmin(2)
                        adt%xBBox(3, nn) = xmin(3)
 
                        adt%xBBox(4, nn) = xmax(1)
                        adt%xBBox(5, nn) = xmax(2)
                        adt%xBBox(6, nn) = xmax(3)
 
                    end if testwithin
                end do
            end do elementloop2
 
            ! Deallocate the memory for elementWithinBBox.
 
            deallocate (elementwithinbbox, stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildVolumeADT", &
                 "Deallocation failure for &
                 &elementWithinBBox.")
 
            else testbbox
 
            ! No global bounding box. The number of local bounding boxes
            ! to be stored is the total number of local volume elements.
 
            ii = ntetra + npyra + nprisms + nhexa
            adt%nBBoxes = ii
 
            ! Allocate the memory for the bounding box coordinates, the
            ! corresponding element type and the index in the connectivity.
 
            allocate (adt%xBBox(6, ii), adt%elementType(ii), &
                      adt%elementID(ii), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "buildVolumeADT", &
                 "Memory allocation failure for bounding &
                 &box data.")
 
            ! Loop over the number of element types present, i.e. 4,
            ! to store the bounding boxes; nn is the counter.
 
            nn = 0
            elementloop3: do ll = 1, 4
 
                ! Set the correct pointers for this element.
 
                call setvolumepointers(ll)
 
                ! Loop over the number of elements and store the bounding
                ! box info.
 
                do i = 1, nelem
                    nn = nn + 1
 
                    adt%elementType(nn) = eltype
                    adt%elementID(nn) = i
 
                    mm = conn(1, i)
                    xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
                    xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
                    xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
                    do j = 2, nnpe
                        mm = conn(j, i)
 
                        xmin(1) = min(xmin(1), coor(1, mm))
                        xmin(2) = min(xmin(2), coor(2, mm))
                        xmin(3) = min(xmin(3), coor(3, mm))
 
                        xmax(1) = max(xmax(1), coor(1, mm))
                        xmax(2) = max(xmax(2), coor(2, mm))
                        xmax(3) = max(xmax(3), coor(3, mm))
                    end do
 
                    adt%xBBox(1, nn) = xmin(1)
                    adt%xBBox(2, nn) = xmin(2)
                    adt%xBBox(3, nn) = xmin(3)
 
                    adt%xBBox(4, nn) = xmax(1)
                    adt%xBBox(5, nn) = xmax(2)
                    adt%xBBox(6, nn) = xmax(3)
                end do
 
            end do elementloop3
 
        end if testbbox
 
        ! Build the ADT from the now known boundary boxes.
 
        call buildadt(adt)
 
        !===============================================================
 
    contains
 
        !=============================================================
 
        subroutine setvolumepointers(ll)
            !
            !          This internal subroutine sets the pointers to the correct
            !          volume element, such that a loop over the element types
            !          can be used.
            !          Subroutine intent(in) arguments.
            !          --------------------------------
            !          ll: Element type for which the pointers must be used.
            !
            implicit none
            !
            !         Subroutine arguments.
            !
            integer, intent(in) :: ll
 
            select case (ll)
            case (1)
                eltype = adttetrahedron; nelem = ntetra; nnpe = 4
                conn => tetraconn
            case (2)
                eltype = adtpyramid; nelem = npyra; nnpe = 5
                conn => pyraconn
            case (3)
                eltype = adtprism; nelem = nprisms; nnpe = 6
                conn => prismsconn
            case (4)
                eltype = adthexahedron; nelem = nhexa; nnpe = 8
                conn => hexaconn
            end select
 
        end subroutine setvolumepointers
 
    end subroutine buildvolumeadt
 
    subroutine buildserialhex(nHexa, nNodes, coor, hexaConn, ADT)
        !
        !        This a specialized routine that builds and ADT tree for
        !        hex volumes only and only in serial. Also, this routine does
        !        use adtDats's ADTs() array list...the user must supply the
        !        adtType to use and is responsible for all data management of
        !        this type.
        !        The memory intensive part of these
        !        arguments, the arrays with the coordinates and
        !        connectivities, are not copied. Instead pointers are set to
        !        these arrays. It is therefore the responsibility of the user
        !        not to deallocate this memory before all the searches have
        !        been performed.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !       * nHexa :    Number of hexa cells
        !        nNodes:    Number of nodes in the given grid.
        !        Subroutine intent(in), target arguments.
        !        ----------------------------------------
        !        coor(3,nNodes):        Nodal coordinates of the local grid.
        !        hexaConn(8,nHexa):     Idem for the hexahedra.
        !        Subroutine intent(out), arguments.
        !        ----------------------------------------
        !        ADT : The newly completed ADT
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nHexa
        integer(kind=intType), intent(in) :: nNodes
 
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: hexaConn
 
        real(kind=realtype), dimension(:, :), intent(in), &
            target :: coor
        type(adttype), intent(out) :: adt
        !
        !       Local variables.
        !
        integer :: ierr, ll, nnpe
        integer(kind=intType) :: i, j, mm
        real(kind=realtype), dimension(3) :: xmin, xmax
 
        ! We need to set comm...explictly mpi_comm_self
        adt%comm = mpi_comm_self
        adt%nProcs = 1
        adt%myID = 0
        ! Set the ADT type, which is a volume ADT.
 
        adt%adtType = adtvolumeadt
 
        ! Copy the number of nodes and volume elements and set the number
        ! of surface elements to 0; only a volume grid has been given.
 
        adt%nNodes = nnodes
        adt%nHexa = nhexa
        adt%nTetra = 0
        adt%nPyra = 0
        adt%nPrisms = 0
        adt%nTria = 0
        adt%nQuads = 0
 
        ! Set the pointers for the coordinates and the
        ! volume connectivities.
 
        adt%coor => coor
        adt%hexaConn => hexaconn
        nullify (adt%tetraConn, adt%pyraConn, adt%prismsConn)
 
        adt%nBBoxes = nhexa
 
        ! Allocate the memory for the bounding box coordinates, the
        ! corresponding element type and the index in the connectivity.
 
        allocate (adt%xBBox(6, nhexa))
        allocate (adt%elementType(nhexa))
        allocate (adt%elementID(nhexa))
 
        ! All hexas
        adt%elementType = adthexahedron
 
        ! Loop over the number of elements and store the bounding
        ! box info.
        nnpe = 8
 
        do i = 1, nhexa
 
            mm = hexaconn(1, i)
            xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
            xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
            xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
            do j = 2, nnpe
                mm = hexaconn(j, i)
 
                xmin(1) = min(xmin(1), coor(1, mm))
                xmin(2) = min(xmin(2), coor(2, mm))
                xmin(3) = min(xmin(3), coor(3, mm))
 
                xmax(1) = max(xmax(1), coor(1, mm))
                xmax(2) = max(xmax(2), coor(2, mm))
                xmax(3) = max(xmax(3), coor(3, mm))
            end do
 
            adt%xBBox(1, i) = xmin(1)
            adt%xBBox(2, i) = xmin(2)
            adt%xBBox(3, i) = xmin(3)
 
            adt%xBBox(4, i) = xmax(1)
            adt%xBBox(5, i) = xmax(2)
            adt%xBBox(6, i) = xmax(3)
 
            ! elementID is just sequential since we only have 1 element type
            adt%elementID(i) = i
 
        end do
 
        ! Build the ADT from the now known boundary boxes.
 
        call buildadt(adt)
 
    end subroutine buildserialhex
 
    subroutine destroyserialhex(ADT)
        ! Deallocate the data allocated from the ADT
 
        implicit none
        type(adttype), intent(inout) :: adt
 
        deallocate (adt%xBBox)
        deallocate (adt%elementType)
        deallocate (adt%elementID)
        deallocate (adt%ADTree)
    end subroutine destroyserialhex
 
    subroutine buildserialquad(nQuad, nNodes, coor, quadsConn, ADT)
        !
        !        This a specialized routine that builds and ADT tree for
        !        quad surface meshes and only in serial. Also, this routine
        !        does not
        !        use adtDats's ADTs() array list...the user must supply the
        !        adtType to use and is responsible for all data management of
        !        this type.
        !        The memory intensive part of these
        !        arguments, the arrays with the coordinates and
        !        connectivities, are not copied. Instead pointers are set to
        !        these arrays. It is therefore the responsibility of the user
        !        not to deallocate this memory before all the searches have
        !        been performed.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !       * nQuad :    Number of quad cells
        !        nNodes:    Number of nodes in the given grid.
        !        Subroutine intent(in), target arguments.
        !        ----------------------------------------
        !        coor(3,nNodes):        Nodal coordinates of the local grid.
        !        quadsConn(8,nHexa):     Connectivity for quad cells
        !        Subroutine intent(out), arguments.
        !        ----------------------------------------
        !        ADT : The newly completed ADT
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nquad
        integer(kind=intType), intent(in) :: nnodes
 
        integer(kind=intType), dimension(:, :), intent(in), &
            target :: quadsconn
 
        real(kind=realtype), dimension(:, :), intent(in), &
            target :: coor
        type(adttype), intent(out) :: adt
        !
        !       Local variables.
        !
        integer :: ierr, ll, nnpe
        integer(kind=intType) :: i, j, mm
        real(kind=realtype), dimension(3) :: xmin, xmax
 
        ! We need to set comm...explictly mpi_comm_self
        adt%comm = mpi_comm_self
        adt%nProcs = 1
        adt%myID = 0
        ! Set the ADT type, which is a surface ADT.
 
        adt%adtType = adtsurfaceadt
 
        ! Copy the number of nodes and volume elements and set the number
        ! of surface elements to 0; only a volume grid has been given.
 
        adt%nNodes = nnodes
        adt%nHexa = 0
        adt%nTetra = 0
        adt%nPyra = 0
        adt%nPrisms = 0
        adt%nTria = 0
        adt%nQuads = nquad
 
        ! Set the pointers for the coordinates and the
        ! volume connectivities.
 
        adt%coor => coor
        adt%quadsConn => quadsconn
        nullify (adt%triaConn)
 
        adt%nBBoxes = nquad
 
        ! Allocate the memory for the bounding box coordinates, the
        ! corresponding element type and the index in the connectivity.
 
        allocate (adt%xBBox(6, nquad))
        allocate (adt%elementType(nquad))
        allocate (adt%elementID(nquad))
 
        ! All hexas
        adt%elementType = adtquadrilateral
 
        ! Loop over the number of elements and store the bounding
        ! box info.
        nnpe = 4
 
        do i = 1, nquad
 
            mm = quadsconn(1, i)
            xmin(1) = coor(1, mm); xmax(1) = coor(1, mm)
            xmin(2) = coor(2, mm); xmax(2) = coor(2, mm)
            xmin(3) = coor(3, mm); xmax(3) = coor(3, mm)
 
            do j = 2, nnpe
                mm = quadsconn(j, i)
 
                xmin(1) = min(xmin(1), coor(1, mm))
                xmin(2) = min(xmin(2), coor(2, mm))
                xmin(3) = min(xmin(3), coor(3, mm))
 
                xmax(1) = max(xmax(1), coor(1, mm))
                xmax(2) = max(xmax(2), coor(2, mm))
                xmax(3) = max(xmax(3), coor(3, mm))
            end do
 
            adt%xBBox(1, i) = xmin(1)
            adt%xBBox(2, i) = xmin(2)
            adt%xBBox(3, i) = xmin(3)
 
            adt%xBBox(4, i) = xmax(1)
            adt%xBBox(5, i) = xmax(2)
            adt%xBBox(6, i) = xmax(3)
 
            ! elementID is just sequential since we only have 1 element type
            adt%elementID(i) = i
 
        end do
 
        ! Build the ADT from the now known boundary boxes.
 
        call buildadt(adt)
 
    end subroutine buildserialquad
 
    subroutine destroyserialquad(ADT)
        ! Deallocate the data allocated from the ADT
 
        implicit none
        type(adttype), intent(inout) :: adt
 
        deallocate (adt%xBBox)
        deallocate (adt%elementType)
        deallocate (adt%elementID)
        deallocate (adt%ADTree)
    end subroutine destroyserialquad
end module adtbuild