surfaceUtils.F90

Go to the documentation of this file.

module surfaceutils
 
contains
 
    subroutine getsurfacesize(sizeNode, sizeCell, famList, n, includeZipper)
        ! Compute the number of points that will be returned from getForces
        ! or getForcePoints
        use constants
        use blockpointers, only: bcdata, ndom, nbocos
        use utils, only: setpointers
        use sorting, only: faminlist
        use surfacefamilies, only: bcfamgroups
        use oversetdata, only: zippermeshes, zippermesh, oversetpresent
 
        implicit none
 
        integer(kind=intType), intent(out) :: sizeNode, sizeCell
        logical, intent(in) :: includeZipper
        integer(kind=intType) :: nn, mm, i, j, iimax, shp(1)
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, iBCGroup
        integer(kind=intType), intent(in) :: n, famList(n)
        type(zippermesh), pointer :: zipper
        logical :: BCGroupNeeded
        sizenode = 0_inttype
        sizecell = 0_inttype
 
        domains: do nn = 1, ndom
            call setpointers(nn, 1_inttype, 1_inttype)
            bocos: do mm = 1, nbocos
                ! Check if this surface should be included or not:
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                    sizenode = sizenode + (iend - ibeg + 1) * (jend - jbeg + 1)
 
                    ! If we don't care about blanking, it's easy:
                    blanking: if (.not. includezipper) then
                        sizecell = sizecell + (iend - ibeg) * (jend - jbeg)
                    else
                        ! Otherwise we have to consider the iBlank
                        do j = jbeg + 1, jend
                            do i = ibeg + 1, iend
                                if (bcdata(mm)%iBlank(i, j) == 1) then
                                    sizecell = sizecell + 1
                                end if
                            end do
                        end do
                    end if blanking
                end if faminclude
            end do bocos
        end do domains
 
        ! We know must consider additional nodes that are required by the
        ! zipper mesh triangles on the root proc.
 
        ! No overset or we don't want to include the zipper, return immediately
        if (.not. oversetpresent .or. .not. includezipper) then
            return
        end if
 
        ! If there are zipper meshes, we must include the nodes that the
        ! zipper triangles will use.
        do ibcgroup = 1, nfamexchange
            bcgroupneeded = .false.
            bcgroupfamloop: do j = 1, size(bcfamgroups(ibcgroup)%famList)
                if (faminlist(bcfamgroups(ibcgroup)%famList(j), famlist)) then
                    bcgroupneeded = .true.
                    exit bcgroupfamloop
                end if
            end do bcgroupfamloop
 
            if (.not. bcgroupneeded) then
                cycle
            end if
 
            ! Pointer for easier reading.
            zipper => zippermeshes(ibcgroup)
 
            ! If we don't have a zipper for this BCGroup, just keep going.
            if (.not. zipper%allocated) then
                cycle
            end if
 
            ! Include the total extra number of nodes. Not necessairly all
            ! nodes are needed, but they will be returned anyway.
            sizenode = sizenode + size(zipper%indices)
 
            ! Include the extra number of cells. Not necessairly all cells
            ! are needed, but here we have to check indvidually.
 
            do i = 1, size(zipper%fam)
                if (faminlist(zipper%fam(i), famlist)) then
                    sizecell = sizecell + 1
                end if
            end do
        end do
 
    end subroutine getsurfacesize
 
    subroutine getsurfaceconnectivity(conn, cgnsBlockID, ncell, famList, nFamList, includeZipper)
        ! Return the connectivity list for the each of the patches
        ! cgnsBlockID is the domain number of the CGNS file that each face patch corresponds to.
        ! the zipper mesh patches will have cgnsBlockID = -1
        use constants
        use blockpointers, only: ndom, nbocos, bcdata, bcfaceid, righthanded, nbkglobal
        use utils, only: setpointers
        use sorting, only: faminlist
        use surfacefamilies, only: bcfamgroups
        use oversetdata, only: zippermeshes, zippermesh, oversetpresent
 
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: ncell
        integer(kind=intType), intent(inout) :: conn(4 * ncell), cgnsBlockID(ncell)
        integer(kind=intType), intent(in) :: nFamList, famList(nFamList)
        logical, intent(in) :: includeZipper
 
        ! Working
        integer(kind=intType) :: nn, mm, cellCount, nodeCount, ni, nj, i, j
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, iBCGroup
        logical regularOrdering, BCGroupNeeded
        type(zippermesh), pointer :: zipper
        cellcount = 0
        nodecount = 0
 
        ! Initialize all cgns IDs to -1. This will take care of the zipper meshes IDs, which should
        ! be -1 since they do not belong to the CGNS file.
        cgnsblockid = -1
 
        domains: do nn = 1, ndom
            call setpointers(nn, 1_inttype, 1_inttype)
            bocos: do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    ni = iend - ibeg + 1
                    nj = jend - jbeg + 1
 
                    ! We want to ensure that all the normals of the faces are
                    ! consistent. To ensure this, we enforce that all normals
                    ! are "into" the domain. Therefore we must treat difference
                    ! faces of a block differently. For example for an iLow
                    ! face, when looping over j-k in the regular way, results
                    ! in in a domain inward pointing normal for iLow but
                    ! outward pointing normal for iHigh. The same is true for
                    ! kMin and kMax. However, it is reverse for the J-faces:
                    ! This is becuase the way the pointers are extracted i then
                    ! k is the reverse of what "should" be for consistency. The
                    ! other two, the pointers are cyclic consistent: i,j->k,
                    ! j,k (wrap) ->i, but for the j-direction is is i,k->j when
                    ! to be consistent with the others it should be
                    ! k,i->j. Hope that made sense.
 
                    select case (bcfaceid(mm))
                    case (imin, jmax, kmin)
                        regularordering = .true.
                    case default
                        regularordering = .false.
                    end select
 
                    ! Now this can be reversed *again* if we have a block that
                    ! is left handed.
                    if (.not. righthanded) then
                        regularordering = .not. (regularordering)
                    end if
 
                    if (regularordering) then
                        ! Do regular ordering.
 
                        ! Loop over generic face size...Note we are doing zero
                        ! based ordering!
 
                        ! This cartoon of a generic cell might help:
                        !
                        ! i, j+1 +-----+ i+1, j+1
                        !   n4   |     | n3
                        !        +-----+
                        !       i,j    i+1, j
                        !       n1     n2
                        !
 
                        do j = 0, nj - 2
                            do i = 0, ni - 2
                                if (.not. includezipper .or. bcdata(mm)%iBlank(i + ibeg + 1, j + jbeg + 1) == 1) then
                                    conn(4 * cellcount + 1) = nodecount + (j) * ni + i + 1! n1
                                    conn(4 * cellcount + 2) = nodecount + (j) * ni + i + 1 + 1! n2
                                    conn(4 * cellcount + 3) = nodecount + (j + 1) * ni + i + 1 + 1! n3
                                    conn(4 * cellcount + 4) = nodecount + (j + 1) * ni + i + 1! n4
 
                                    ! Assign the corresponding block IDs.
                                    ! Remember that nbkGlobal is the CGNS ID of the current domain,
                                    ! and this is set when we call setPointers.
                                    cgnsblockid(cellcount + 1) = nbkglobal
 
                                    cellcount = cellcount + 1
                                end if
                            end do
                        end do
                    else
                        ! Do reverse ordering:
                        do j = 0, nj - 2
                            do i = 0, ni - 2
                                if (.not. includezipper .or. bcdata(mm)%iBlank(i + ibeg + 1, j + jbeg + 1) == 1) then
                                    conn(4 * cellcount + 1) = nodecount + (j) * ni + i + 1! n1
                                    conn(4 * cellcount + 2) = nodecount + (j + 1) * ni + i + 1! n4
                                    conn(4 * cellcount + 3) = nodecount + (j + 1) * ni + i + 1 + 1! n3
                                    conn(4 * cellcount + 4) = nodecount + (j) * ni + i + 1 + 1! n2
 
                                    ! Assign the corresponding block IDs.
                                    ! Remember that nbkGlobal is the CGNS ID of the current domain,
                                    ! and this is set when we call setPointers.
                                    cgnsblockid(cellcount + 1) = nbkglobal
 
                                    cellcount = cellcount + 1
                                end if
                            end do
                        end do
                    end if
                    nodecount = nodecount + ni * nj
                end if faminclude
            end do bocos
        end do domains
 
        ! We know must consider additional connectivity required by the
        ! zipper mesh triangles on the root proc
 
        ! No overset or don't want zipper return immediately
        if (.not. oversetpresent .or. .not. includezipper) then
            return
        end if
 
        ! If there are zipper meshes, we must include the nodes that the
        ! zipper triangles will use.
        bcgrouploop: do ibcgroup = 1, nfamexchange
 
            bcgroupneeded = .false.
            bcgroupfamloop: do i = 1, size(bcfamgroups(ibcgroup)%famList)
                if (faminlist(bcfamgroups(ibcgroup)%famList(i), famlist)) then
                    bcgroupneeded = .true.
                    exit bcgroupfamloop
                end if
            end do bcgroupfamloop
 
            if (.not. bcgroupneeded) then
                cycle
            end if
 
            ! Pointer for easier reading.
            zipper => zippermeshes(ibcgroup)
 
            ! If the zipper isn't done yet, don't do anything
            if (.not. zipper%allocated) then
                cycle
            end if
 
            ! Include the extra number of cells. Not necessairly all cells
            ! are needed, but ehre we have to check indvidually.
 
            do i = 1, size(zipper%fam)
                if (faminlist(zipper%fam(i), famlist)) then
                    ! This triangle should be included. Note that we use
                    ! degenerate quads for the triangles.
                    conn(4 * cellcount + 1) = nodecount + zipper%conn(1, i)
                    conn(4 * cellcount + 2) = nodecount + zipper%conn(2, i)
                    conn(4 * cellcount + 3) = nodecount + zipper%conn(3, i)
                    conn(4 * cellcount + 4) = nodecount + zipper%conn(3, i)
                    cellcount = cellcount + 1
                end if
            end do
        end do bcgrouploop
 
    end subroutine getsurfaceconnectivity
 
    subroutine getsurfacefamily(elemFam, ncell, famList, nFamList, includeZipper)
 
        use constants
        use blockpointers, only: ndom, nbocos, bcdata
        use utils, only: setpointers
        use sorting, only: faminlist
        use surfacefamilies, only: bcfamgroups
        use oversetdata, only: zippermeshes, zippermesh, oversetpresent
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: ncell
        integer(kind=intType), intent(inout) :: elemFam(nCell)
        integer(kind=intType), intent(in) :: nFamList, famList(nFamList)
        logical, intent(in) :: includeZipper
 
        ! Working
        integer(kind=intType) :: nn, mm, cellCount, nodeCount, ni, nj, i, j
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, iBCGroup
        logical BCGroupNeeded
        type(zippermesh), pointer :: zipper
 
        cellcount = 0
        nodecount = 0
 
        domains: do nn = 1, ndom
            call setpointers(nn, 1_inttype, 1_inttype)
            bocos: do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    ni = iend - ibeg + 1
                    nj = jend - jbeg + 1
                    do j = 0, nj - 2
                        do i = 0, ni - 2
                            if (.not. includezipper .or. bcdata(mm)%iBlank(i + ibeg + 1, j + jbeg + 1) == 1) then
                                cellcount = cellcount + 1
                                elemfam(cellcount) = bcdata(mm)%famID
                            end if
                        end do
                    end do
                    nodecount = nodecount + ni * nj
                end if faminclude
            end do bocos
        end do domains
 
        ! We know must consider additional elements quired by the zipper
        ! mesh triangles on the root proc
 
        ! No overset or don't want zipper
        if (.not. oversetpresent .or. .not. includezipper) then
            return
        end if
 
        ! If there are zipper meshes, we must include the nodes that the
        ! zipper triangles will use.
        bcgrouploop: do ibcgroup = 1, nfamexchange
 
            bcgroupneeded = .false.
            bcgroupfamloop: do i = 1, size(bcfamgroups(ibcgroup)%famList)
                if (faminlist(bcfamgroups(ibcgroup)%famList(i), famlist)) then
                    bcgroupneeded = .true.
                    exit bcgroupfamloop
                end if
            end do bcgroupfamloop
 
            if (.not. bcgroupneeded) then
                cycle
            end if
 
            ! Pointer for easier reading.
            zipper => zippermeshes(ibcgroup)
 
            ! If the zipper isn't done yet, don't do anything
            if (.not. zipper%allocated) then
                cycle
            end if
 
            ! Include the extra number of cells. Not necessairly all cells
            ! are needed, but ehre we have to check indvidually.
 
            do i = 1, size(zipper%fam)
                if (faminlist(zipper%fam(i), famlist)) then
                    ! This triangle should be included. Note that we use
                    ! degenerate quads for the triangles.
                    cellcount = cellcount + 1
                    elemfam(cellcount) = zipper%fam(i)
                end if
            end do
        end do bcgrouploop
    end subroutine getsurfacefamily
 
    subroutine getsurfacepoints(points, npts, sps_in, famList, nFamList, includeZipper)
        use constants
        use communication, only: myid
        use blockpointers, only: ndom, bcdata, nbocos, x, bcfaceid, il, jl, kl
        use bcpointers, only: xx
        use surfacefamilies, only: bcfamgroups, familyexchange, bcfamexchange
        use oversetdata, only: zippermeshes, zippermesh, oversetpresent
        use sorting, only: faminlist
        use utils, only: setpointers, echk, setbcpointers
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        !
        !      Local variables.
        !
        integer(kind=intType), intent(in) :: npts, sps_in
        real(kind=realtype), intent(inout) :: points(3, npts)
        integer(kind=intType), intent(in) :: nFamList, famList(nFamList)
        logical, intent(in) :: includeZipper
 
        integer(kind=intType) :: mm, nn, i, j, ii, sps, iDim, jj, ierr
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, iBCGroup
        type(zippermesh), pointer :: zipper
        type(familyexchange), pointer :: exch
        logical :: BCGroupNeeded
        real(kind=realtype), dimension(:), pointer :: localptr
        sps = sps_in
 
        ii = 0
        domains: do nn = 1, ndom
            call setpointers(nn, 1_inttype, sps)
 
            ! Loop over the number of boundary subfaces of this block.
            bocos: do mm = 1, nbocos
 
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    ! NODE Based
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    do j = jbeg, jend ! This is a node loop
                        do i = ibeg, iend ! This is a node loop
                            ii = ii + 1
                            select case (bcfaceid(mm))
 
                            case (imin)
                                points(:, ii) = x(1, i, j, :)
                            case (imax)
                                points(:, ii) = x(il, i, j, :)
                            case (jmin)
                                points(:, ii) = x(i, 1, j, :)
                            case (jmax)
                                points(:, ii) = x(i, jl, j, :)
                            case (kmin)
                                points(:, ii) = x(i, j, 1, :)
                            case (kmax)
                                points(:, ii) = x(i, j, kl, :)
                            end select
                        end do
                    end do
                end if faminclude
            end do bocos
        end do domains
 
        ! No overset or not zipper, return
        if (.not. oversetpresent .or. .not. includezipper) then
            return
        end if
 
        ! If there are zipper meshes, we must include the nodes that the
        ! zipper triangles will use.
        do ibcgroup = 1, nfamexchange
 
            zipper => zippermeshes(ibcgroup)
 
            if (.not. zipper%allocated) then
                cycle
            end if
 
            exch => bcfamexchange(ibcgroup, sps)
            bcgroupneeded = .false.
            bcgroupfamloop: do i = 1, size(bcfamgroups(ibcgroup)%famList)
                if (faminlist(bcfamgroups(ibcgroup)%famList(i), famlist)) then
                    bcgroupneeded = .true.
                    exit bcgroupfamloop
                end if
            end do bcgroupfamloop
 
            if (.not. bcgroupneeded) then
                cycle
            end if
 
            ! Now we know we *actually* need something from this BCGroup.
 
            ! Loop over each dimension individually since we have a scalar
            ! scatter.
            dimloop: do idim = 1, 3
 
                call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
                call echk(ierr, __file__, __line__)
                localptr = zero
 
                ! jj is the running counter through the pointer array.
                jj = 0
                do nn = 1, ndom
                    call setpointers(nn, 1_inttype, sps)
                    do mm = 1, nbocos
                        faminclude2: if (faminlist(bcdata(mm)%famID, exch%famList)) then
                            ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                            jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                            call setbcpointers(mm, .true.)
                            do j = jbeg, jend
                                do i = ibeg, iend
                                    jj = jj + 1
                                    localptr(jj) = xx(i + 1, j + 1, idim)
                                end do
                            end do
                        end if faminclude2
                    end do
                end do
 
                ! Restore the pointer
                call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Now scatter this to the zipper
                call vecscatterbegin(zipper%scatter, exch%nodeValLocal, &
                                     zipper%localVal, insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                call vecscatterend(zipper%scatter, exch%nodeValLocal, &
                                   zipper%localVal, insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                ! The values we need are precisely what is in zipper%localVal
                call vecgetarrayf90(zipper%localVal, localptr, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Just copy the received data into the points array. Only root proc.
                if (myid == 0) then
                    points(idim, ii + 1:ii + size(localptr)) = localptr
                end if
                ! The values we need are precisely what is in zipper%localVal
                call vecrestorearrayf90(zipper%localVal, localptr, ierr)
                call echk(ierr, __file__, __line__)
 
            end do dimloop
 
            ! Increcment the running ii counter.
            ii = ii + size(localptr)
 
        end do
    end subroutine getsurfacepoints
 
    subroutine mapvector(vec1, n1, famList1, nf1, vec2, n2, famList2, nf2, includeZipper, nDim)
 
        ! Map one vector, vec1 of size (3,n1) defined on family list 'famList1' onto
        ! vector, vec2, of size (3, n2) defined on family list 'famList2'
 
        ! This operation is actually pretty fast since it just requires a
        ! single copy of surface-based data.
        use constants
        use blockpointers, onlY: ndom, flowdoms
        use sorting, only: faminlist
        use surfacefamilies, only: bcfamgroups
        use oversetdata, only: zippermeshes, zippermesh, oversetpresent
 
        implicit none
 
        ! Input/Output
        integer(kind=intType) :: n1, n2, nf1, nf2, nDim
        integer(kind=intType), intent(in) :: famList1(nf1), famList2(nf2)
        real(kind=realtype), intent(in) :: vec1(ndim, n1)
        real(kind=realtype), intent(inout) :: vec2(ndim, n2)
        logical, intent(in) :: includeZipper
 
        ! Working
        integer(kind=intType) :: i, k, ii, jj, nn, mm, iSize
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, famID, iBCGroup
        logical :: fam1Included, fam2Included
        type(zippermesh), pointer :: zipper
        logical :: BCGroupNeeed
        ii = 0
        jj = 0
        domains: do nn = 1, ndom
            ! Don't set pointers for speed
 
            ! Loop over the number of boundary subfaces of this block.
            bocos: do mm = 1, flowdoms(nn, 1, 1)%nBocos
                famid = flowdoms(nn, 1, 1)%BCdata(mm)%famID
 
                fam1included = faminlist(famid, famlist1)
                fam2included = faminlist(famid, famlist2)
 
                jbeg = flowdoms(nn, 1, 1)%bcData(mm)%jnBeg
                jend = flowdoms(nn, 1, 1)%bcData(mm)%jnEnd
 
                ibeg = flowdoms(nn, 1, 1)%bcData(mm)%inBeg
                iend = flowdoms(nn, 1, 1)%bcData(mm)%inEnd
                isize = (iend - ibeg + 1) * (jend - jbeg + 1)
 
                if (fam1included .and. fam2included) then
                    ! The two lists overlap so copy:
                    do k = 1, isize
                        vec2(:, k + jj) = vec1(:, k + ii)
                    end do
                end if
 
                ! Finally increment the counters if the face had been inclded
                if (fam1included) then
                    ii = ii + isize
                end if
 
                if (fam2included) then
                    jj = jj + isize
                end if
 
            end do bocos
        end do domains
 
        ! As with the rest of the code we have to account for the zipper
        ! mesh on the root proc.
 
        ! We know must consider additional nodes that are required by the
        ! zipper mesh triangles on the root proc.
 
        ! No overset or don't want to include zipper, return immediately
        if (.not. oversetpresent .or. .not. includezipper) then
            return
        end if
 
        ! If there are zipper meshes, we must include the nodes that the
        ! zipper triangles will use.
        bcgrouploop: do ibcgroup = 1, nfamexchange
 
            zipper => zippermeshes(ibcgroup)
 
            fam1included = .false.
            fam2included = .false.
            bcgroupfamloop: do i = 1, size(bcfamgroups(ibcgroup)%famList)
                if (faminlist(bcfamgroups(ibcgroup)%famList(i), famlist1)) then
                    fam1included = .true.
                end if
                if (faminlist(bcfamgroups(ibcgroup)%famList(i), famlist2)) then
                    fam2included = .true.
                end if
            end do bcgroupfamloop
 
            ! This is the total number of nodes that this BCGroup has. It
            ! is not further broken down by family group.
            isize = size(zipper%indices)
 
            if (fam1included .and. fam2included) then
                ! The two lists overlap so copy:
                do k = 1, isize
                    if (k + ii <= n1) then
                        vec2(:, k + jj) = vec1(:, k + ii)
                    end if
                end do
            end if
 
            ! Finally increment the counters if this BCGroup had been included.
            if (fam1included) then
                ii = ii + isize
            end if
 
            if (fam2included) then
                jj = jj + isize
            end if
        end do bcgrouploop
 
    end subroutine mapvector
 
    subroutine getwalllist(wallList, nWallList, nFamTotal)
 
        ! Python wrapped utility function to return the list of families
        ! that are walls to Python since we need that information in
        ! Python for a few default values.
 
        use constants
        use surfacefamilies, only: bcfamgroups
        implicit none
 
        integer(kind=intType), intent(in) :: nFamtotal
        integer(kind=intType), dimension(nFamTotal), intent(out) :: wallList
        integer(kind=intType), intent(out) :: nWallList
 
        nwalllist = size(bcfamgroups(ibcgroupwalls)%famList)
        walllist(1:nwalllist) = bcfamgroups(ibcgroupwalls)%famList
 
    end subroutine getwalllist
 
end module surfaceutils