getForces.F90

Go to the documentation of this file.

 
subroutine getforces(forces, npts, sps)
    use constants
    use communication, only: myid
    use blockpointers, only: bcdata, ndom, nbocos, bctype
    use inputphysics, only: forcesastractions
    use utils, only: setpointers, terminate, echk
    use surfaceintegrations, only: integratesurfaces
    use surfacefamilies, only: fullfamlist
    use oversetdata, only: zippermeshes, zippermesh, oversetpresent
    use surfacefamilies, only: familyexchange, bcfamexchange
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
    integer(kind=intType), intent(in) :: npts, sps
    real(kind=realtype), intent(inout) :: forces(3, npts)
 
    integer(kind=intType) :: mm, nn, i, j, ii, jj, iDim, ierr
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
    real(kind=realtype) :: sss(3), v2(3), v1(3), qa, sepsensor, cavitation
    real(kind=realtype) :: sepsensoravg(3)
    real(kind=realtype) :: fp(3), fv(3), mp(3), mv(3), yplusmax, qf(3)
    real(kind=realtype) :: localvalues(nlocalvalues)
    type(zippermesh), pointer :: zipper
    type(familyexchange), pointer :: exch
    real(kind=realtype), dimension(:), pointer :: localptr
    real(kind=realtype), dimension(nCostFunction) :: funcvalues
    ! Make sure *all* forces are computed. Sectioning will be done
    ! else-where.
 
    domains: do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        localvalues = zero
        call integratesurfaces(localvalues, fullfamlist)
    end do domains
 
    if (forcesastractions) then
        ! Compute tractions if necessary
        call computenodaltractions(sps)
    else
        call computenodalforces(sps)
    end if
 
    ii = 0
    domains2: 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
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
 
                ! This is easy, just copy out F or T in continuous ordering.
                do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd
                    do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd
                        ii = ii + 1
                        if (forcesastractions) then
                            forces(:, ii) = bcdata(mm)%Tp(i, j, :) + bcdata(mm)%Tv(i, j, :)
                        else
                            forces(:, ii) = bcdata(mm)%F(i, j, :)
                        end if
                    end do
                end do
            end if
        end do bocos
    end do domains2
 
    ! We know must consider additional forces that are required by the
    ! zipper mesh triangles on the root proc.
 
    ! Pointer for easier reading.
    zipper => zippermeshes(ibcgroupwalls)
    exch => bcfamexchange(ibcgroupwalls, sps)
    ! No overset present or the zipper isn't allocated nothing to do:
    if (.not. oversetpresent .or. .not. zipper%allocated) then
        return
    end if
 
    if (.not. forcesastractions) then
        ! We have a zipper and regular forces are requested. This is not yet supported.
        call terminate('getForces', 'getForces() is not implmented for zipper meshes and '&
             &'forcesAsTractions=False')
    end if
 
    ! 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__)
 
        ! Copy in the values we already have to the exchange.
        ii = size(localptr)
        localptr = forces(idim, 1:ii)
 
        ! 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 forces array. Just on root proc:
        if (myid == 0) then
            forces(idim, ii + 1:ii + size(localptr)) = localptr
        end if
 
        call vecgetarrayf90(zipper%localVal, localptr, ierr)
        call echk(ierr, __file__, __line__)
    end do dimloop
 
end subroutine getforces
 
subroutine getforces_d(forces, forcesd, npts, sps)
 
    ! This routine performs the forward mode linearization getForces. It
    ! takes in perturbations defined on bcData(mm)%Fp, bcData(mm)%Fv and
    ! bcData(mm)%area and computes either the nodal forces or nodal
    ! tractions.
    use constants
    use communication, only: myid
    use blockpointers, only: ndom, nbocos, bcdata, bctype, nbocos, bcdatad
    use inputphysics, only: forcesastractions
    use surfacefamilies, only: bcfamexchange, familyexchange
    use utils, only: setpointers, setpointers_d, echk, terminate
    use oversetdata, only: zippermeshes, zippermesh, oversetpresent
    use surfacefamilies, only: familyexchange, bcfamexchange
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
    integer(kind=intType), intent(in) :: npts, sps
    real(kind=realtype), intent(out), dimension(3, npts) :: forces, forcesd
    integer(kind=intType) :: mm, nn, i, j, ii, jj, iDim, ierr
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, ind(4), ni, nj
    real(kind=realtype) :: qa, qad, qf, qfd
    real(kind=realtype), dimension(:), pointer :: localptr, localptrd
    type(zippermesh), pointer :: zipper
    type(familyexchange), pointer :: exch
 
    if (forcesastractions) then
        call computenodaltractions_d(sps)
    else
        call computenodalforces_d(sps)
    end if
 
    ! Extract the values out into the output derivative array
    ii = 0
    domains2: do nn = 1, ndom
        call setpointers_d(nn, 1_inttype, sps)
 
        ! Loop over the number of boundary subfaces of this block.
        bocos: do mm = 1, nbocos
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
 
                ! This is easy, just copy out F or T in continuous ordering.
                do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd
                    do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd
                        ii = ii + 1
                        if (forcesastractions) then
                            forcesd(:, ii) = bcdatad(mm)%Tp(i, j, :) + bcdatad(mm)%Tv(i, j, :)
                        else
                            forcesd(:, ii) = bcdatad(mm)%F(i, j, :)
                        end if
                    end do
                end do
            end if
        end do bocos
    end do domains2
 
    ! We know must consider additional forces that are required by the
    ! zipper mesh triangles on the root proc.
 
    ! Pointer for easier reading.
    zipper => zippermeshes(ibcgroupwalls)
    exch => bcfamexchange(ibcgroupwalls, sps)
    ! No overset present or the zipper isn't allocated nothing to do:
    if (.not. oversetpresent .or. .not. zipper%allocated) then
        return
    end if
 
    if (.not. forcesastractions) then
        ! We have a zipper and regular forces are requested. This is not yet supported.
        call terminate('getForces', 'getForces() is not implmented for zipper meshes and '&
             &'forcesAsTractions=False')
    end if
 
    ! 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__)
 
        ! Copy in the values we already have to the exchange.
        ii = size(localptr)
        localptr = forcesd(idim, 1:ii)
 
        ! 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 forces array. Just on root proc:
        if (myid == 0) then
            forcesd(idim, ii + 1:ii + size(localptr)) = localptr
        end if
 
        call vecgetarrayf90(zipper%localVal, localptr, ierr)
        call echk(ierr, __file__, __line__)
    end do dimloop
 
end subroutine getforces_d
 
subroutine getforces_b(forcesd, npts, sps)
 
    ! This routine performs the reverse of getForces. It takes in
    ! forces_b and perfroms the reverse of the nodal averaging procedure
    ! in getForces to compute bcDatad(mm)%Fp, bcDatad(mm)%Fv and
    ! bcDatad(mm)%area.
    use constants
    use communication, only: myid
    use blockpointers, only: ndom, nbocos, bcdata, bctype, nbocos, bcdatad
    use inputphysics, only: forcesastractions
    use surfacefamilies, only: bcfamexchange, familyexchange
    use utils, only: echk, setpointers, setpointers_d
    use oversetdata, only: zippermeshes, zippermesh, oversetpresent
    use surfacefamilies, only: familyexchange, bcfamexchange
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
    integer(kind=intType), intent(in) :: npts, sps
    real(kind=realtype), intent(inout) :: forcesd(3, npts)
    integer(kind=intType) :: mm, nn, i, j, ii, iDim, ierr
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
    type(zippermesh), pointer :: zipper
    type(familyexchange), pointer :: exch
    real(kind=realtype), dimension(:), pointer :: localptr
    real(kind=realtype), dimension(3, npts) :: forces
 
    ! Run nonlinear code to make sure that all intermediate values are
    ! updated.
    call getforces(forces, npts, sps)
 
    ! We know must consider additional forces that are required by the
    ! zipper mesh triangles on the root proc.
 
    ! Pointer for easier reading.
    zipper => zippermeshes(ibcgroupwalls)
    exch => bcfamexchange(ibcgroupwalls, sps)
    ! No overset present or the zipper isn't allocated nothing to do:
    zipperreverse: if (oversetpresent .and. zipper%allocated) then
 
        ! Loop over each dimension individually since we have a scalar
        ! scatter.
        dimloop: do idim = 1, 3
 
            call vecgetarrayf90(zipper%localVal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            ii = exch%nNodes
            ! Just copy the received data into the forces array. Just on root proc:
            if (myid == 0) then
                do i = 1, size(localptr)
                    localptr(i) = forcesd(idim, ii + i)
                    forcesd(idim, ii + i) = zero
                end do
            end if
 
            call vecgetarrayf90(zipper%localVal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Zero the vector we are scatting into:
            call vecset(exch%nodeValLocal, zero, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Scatter values from the root using the zipper scatter.
            call vecscatterbegin(zipper%scatter, zipper%localVal, &
                                 exch%nodeValLocal, add_values, scatter_reverse, ierr)
            call echk(ierr, __file__, __line__)
 
            call vecscatterend(zipper%scatter, zipper%localVal, &
                               exch%nodeValLocal, add_values, scatter_reverse, ierr)
            call echk(ierr, __file__, __line__)
 
            call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Accumulate the scatted values onto forcesd
            ii = size(localptr)
            forcesd(idim, 1:ii) = forcesd(idim, 1:ii) + localptr
 
            ! Restore the pointer
            call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
        end do dimloop
    end if zipperreverse
 
    ! Set the incoming derivative values
    ii = 0
    domains2: do nn = 1, ndom
        call setpointers_d(nn, 1_inttype, sps)
 
        ! Loop over the number of boundary subfaces of this block.
        bocos: do mm = 1, nbocos
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
                ! This is easy, just copy out F or T in continuous ordering.
                do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd
                    do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd
                        ii = ii + 1
                        if (forcesastractions) then
                            bcdatad(mm)%Tp(i, j, :) = forcesd(:, ii)
                            bcdatad(mm)%Tv(i, j, :) = forcesd(:, ii)
                        else
                            bcdatad(mm)%F(i, j, :) = forcesd(:, ii)
                        end if
                    end do
                end do
            end if
        end do bocos
    end do domains2
 
    if (.not. forcesastractions) then
        ! For forces, we can accumulate the nodal seeds on the Fp and Fv
        ! values. The area seed is zeroed.
        call computenodalforces_b(sps)
    else
        call computenodaltractions_b(sps)
    end if
 
end subroutine getforces_b
 
subroutine surfacecellcentertonode(exch)
 
    use constants
    use blockpointers, only: bcdata, ndom, nbocos, bctype
    use surfacefamilies, only: familyexchange
    use utils, only: setpointers, echk
    use sorting, only: faminlist
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
 
    type(familyexchange) :: exch
    integer(kind=intType) :: sps
    integer(kind=intType) :: mm, nn, i, j, ii, jj, idim, ierr
    integer(kind=intType) :: ibeg, iend, jbeg, jend, ind(4), ni, nj
    real(kind=realtype) :: qv
    real(kind=realtype), dimension(:), pointer :: localptr
 
    ! We assume that normalization factor is already computed
    sps = exch%sps
    call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
    localptr = zero
 
    ! ii is the running counter through the pointer array.
    ii = 0
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            faminclude: if (faminlist(bcdata(mm)%famID, exch%famList)) then
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                do j = 0, nj - 2
                    do i = 0, ni - 2
                        ! Note: No +iBeg, and +jBeg becuase cellVal is a pointer
                        ! and always starts at one
                        qv = fourth * bcdata(mm)%cellVal(i + 1, j + 1)
                        ind(1) = ii + (j) * ni + i + 1
                        ind(2) = ii + (j) * ni + i + 2
                        ind(3) = ii + (j + 1) * ni + i + 2
                        ind(4) = ii + (j + 1) * ni + i + 1
                        do jj = 1, 4
                            localptr(ind(jj)) = localptr(ind(jj)) + qv
                        end do
                    end do
                end do
                ii = ii + ni * nj
            end if faminclude
        end do
    end do
 
    call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Globalize the current face based value
    call vecset(exch%nodeValGlobal, zero, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                         exch%nodeValGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, exch%nodeValLocal, &
                       exch%nodeValGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Now divide by the weighting. We can do this with a vecpointwisemult
    call vecpointwisemult(exch%nodeValGlobal, exch%nodeValGlobal, &
                          exch%sumGlobal, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Push back to the local values
    call vecscatterbegin(exch%scatter, exch%nodeValGlobal, &
                         exch%nodeValLocal, insert_values, scatter_reverse, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, exch%nodeValGlobal, &
                       exch%nodeValLocal, insert_values, scatter_reverse, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ii = 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
 
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                do j = 1, nj
                    do i = 1, ni
                        ! Note: No +iBeg, and +jBeg becuase cellVal is a pointer
                        ! and always starts at one
                        ii = ii + 1
                        bcdata(mm)%nodeVal(i, j) = localptr(ii)
                    end do
                end do
            end if faminclude2
        end do
    end do
 
    call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
end subroutine surfacecellcentertonode
 
subroutine computeweighting(exch)
 
    use constants
    use blockpointers, only: bcdata, ndom, nbocos, bctype
    use surfacefamilies, only: familyexchange
    use utils, only: setpointers, echk
    use sorting, only: faminlist
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
    type(familyexchange) :: exch
    integer(kind=intType) :: sps
    integer(kind=intType) :: mm, nn, i, j, ii, jj, idim, ierr
    integer(kind=intType) :: ibeg, iend, jbeg, jend, ind(4), ni, nj
    real(kind=realtype) :: qf, qa
    real(kind=realtype), dimension(:), pointer :: localptr
 
    sps = exch%sps
 
    call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    localptr = zero
    ! ii is the running counter through the pointer array.
    ii = 0
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            faminclude: if (faminlist(bcdata(mm)%famID, exch%famList)) then
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                do j = 0, nj - 2
                    do i = 0, ni - 2
 
                        ! Scatter a quarter of the face value to each node:
                        ! Note: No +iBeg, and +jBeg becuase cellVal is a pointer
                        ! and always starts at one
                        qa = fourth * bcdata(mm)%cellVal(i + 1, j + 1)
                        ind(1) = ii + (j) * ni + i + 1
                        ind(2) = ii + (j) * ni + i + 2
                        ind(3) = ii + (j + 1) * ni + i + 2
                        ind(4) = ii + (j + 1) * ni + i + 1
                        do jj = 1, 4
                            localptr(ind(jj)) = localptr(ind(jj)) + qa
                        end do
                    end do
                end do
                ii = ii + ni * nj
            end if faminclude
        end do
    end do
 
    call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Globalize the face value
    call vecset(exch%sumGlobal, zero, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                         exch%sumGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, exch%nodeValLocal, &
                       exch%sumGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Now compute the inverse of the weighting so that we can multiply
    ! instead of dividing. Note that we check dividing by zero and just
    ! set those to zero.
 
    call vecgetarrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
    do i = 1, size(localptr)
        if (localptr(i) == zero) then
            localptr(i) = zero
        else
            localptr(i) = one / localptr(i)
        end if
    end do
 
    call vecrestorearrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
end subroutine computeweighting
 
subroutine computenodaltractions(sps)
    use constants
    use blockpointers, only: bcdata, ndom, nbocos, bctype
    use surfacefamilies, only: bcfamexchange, familyexchange
    use utils, only: setpointers, echk, iswalltype
    implicit none
 
    integer(kind=intType), intent(in) :: sps
    integer(kind=intType) :: mm, nn, i, j, ii, jj, iDim, ierr
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, ind(4), ni, nj
    real(kind=realtype) :: qf, qa
    real(kind=realtype), dimension(:), pointer :: localptr
    type(familyexchange), pointer :: exch
 
    ! Set the pointer to the wall exchange:
    exch => bcfamexchange(ibcgroupwalls, sps)
 
    ! Set the weighting factors. In this case, area
    ii = 0
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
 
            bocotype1: if (iswalltype(bctype(mm))) then
                bcdata(mm)%cellVal => bcdata(mm)%area(:, :)
            end if bocotype1
        end do
    end do
    call computeweighting(exch)
 
    fpfvloop: do idim = 1, 6
        ! ii is the running counter through the pointer array.
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, 1_inttype, sps)
            do mm = 1, nbocos
                bocotype2: if (iswalltype(bctype(mm))) then
                    if (idim <= 3) then
                        bcdata(mm)%cellVal => bcdata(mm)%Fp(:, :, idim)
                        bcdata(mm)%nodeVal => bcdata(mm)%Tp(:, :, idim)
                    else
                        bcdata(mm)%cellVal => bcdata(mm)%Fv(:, :, idim - 3)
                        bcdata(mm)%nodeVal => bcdata(mm)%Tv(:, :, idim - 3)
                    end if
                end if bocotype2
            end do
        end do
 
        call surfacecellcentertonode(exch)
 
    end do fpfvloop
 
end subroutine computenodaltractions
 
subroutine computenodaltractions_d(sps)
 
    ! Forward mode lineariation of nodal tractions
 
    use constants
    use blockpointers, only: ndom, nbocos, bcdata, bctype, nbocos, bcdatad
    use inputphysics, only: forcesastractions
    use surfacefamilies, only: bcfamexchange, familyexchange
    use utils, only: setpointers, setpointers_d, echk
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
    integer(kind=intType), intent(in) :: sps
    integer(kind=intType) :: mm, nn, i, j, ii, jj, idim, ierr
    integer(kind=intType) :: ibeg, iend, jbeg, jend, ind(4), ni, nj
    real(kind=realtype) :: qa, qad, qf, qfd
    real(kind=realtype), dimension(:), pointer :: localptr, localptrd
    type(familyexchange), pointer :: exch
    vec nodevallocald, nodevalglobald, sumglobald, tmp
 
    exch => bcfamexchange(ibcgroupwalls, sps)
 
    call vecduplicate(exch%nodeValLocal, nodevallocald, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(exch%nodeValGlobal, nodevalglobald, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(exch%sumGlobal, sumglobald, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(exch%sumGlobal, tmp, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecgetarrayf90(nodevallocald, localptrd, ierr)
    call echk(ierr, __file__, __line__)
 
    localptrd = zero
    localptr = zero
    ! ii is the running counter through the pointer array.
    ii = 0
    do nn = 1, ndom
        call setpointers_d(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
            ni = iend - ibeg + 1
            nj = jend - jbeg + 1
 
            if (bctype(mm) == eulerwall .or. &
                bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
 
                do j = 0, nj - 2
                    do i = 0, ni - 2
 
                        ! Scatter a quarter of the area to each node:
                        qa = fourth * bcdata(mm)%area(i + ibeg + 1, j + jbeg + 1)
                        qad = fourth * bcdatad(mm)%area(i + ibeg + 1, j + jbeg + 1)
                        ind(1) = ii + (j) * ni + i + 1
                        ind(2) = ii + (j) * ni + i + 2
                        ind(3) = ii + (j + 1) * ni + i + 2
                        ind(4) = ii + (j + 1) * ni + i + 1
                        do jj = 1, 4
                            localptrd(ind(jj)) = localptrd(ind(jj)) + qad
                            localptr(ind(jj)) = localptr(ind(jj)) + qa
                        end do
                    end do
                end do
                ii = ii + ni * nj
            end if
        end do
    end do
    call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecrestorearrayf90(nodevallocald, localptrd, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Globalize the area
    call vecset(exch%sumGlobal, zero, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                         exch%sumGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, exch%nodeValLocal, &
                       exch%sumGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Globalize the area derivative
    call vecset(sumglobald, zero, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterbegin(exch%scatter, nodevallocald, &
                         sumglobald, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, nodevallocald, &
                       sumglobald, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Now compute the inverse of the weighting so that we can multiply
    ! instead of dividing. Here we need the original value too:
 
    call vecgetarrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecgetarrayf90(sumglobald, localptrd, ierr)
    call echk(ierr, __file__, __line__)
 
    localptrd = -(localptrd / localptr**2)
    localptr = one / localptr
 
    call vecgetarrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecrestorearrayf90(sumglobald, localptrd, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Now do each of the three dimensions for the pressure and viscous forces
    dimloop: do idim = 1, 6
 
        call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecgetarrayf90(nodevallocald, localptrd, ierr)
        call echk(ierr, __file__, __line__)
 
        localptr = zero
        localptrd = zero
 
        ! ii is the running counter through the pointer array.
        ii = 0
        do nn = 1, ndom
            call setpointers_d(nn, 1_inttype, sps)
            do mm = 1, nbocos
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                if (bctype(mm) == eulerwall .or. &
                    bctype(mm) == nswalladiabatic .or. &
                    bctype(mm) == nswallisothermal) then
                    do j = 0, nj - 2
                        do i = 0, ni - 2
                            if (idim <= 3) then
                                qf = fourth * bcdata(mm)%Fp(i + ibeg + 1, j + jbeg + 1, idim)
                                qfd = fourth * bcdatad(mm)%Fp(i + ibeg + 1, j + jbeg + 1, idim)
                            else
                                qf = fourth * bcdata(mm)%Fv(i + ibeg + 1, j + jbeg + 1, idim - 3)
                                qfd = fourth * bcdatad(mm)%Fv(i + ibeg + 1, j + jbeg + 1, idim - 3)
                            end if
 
                            ind(1) = ii + (j) * ni + i + 1
                            ind(2) = ii + (j) * ni + i + 2
                            ind(3) = ii + (j + 1) * ni + i + 2
                            ind(4) = ii + (j + 1) * ni + i + 1
                            do jj = 1, 4
                                localptr(ind(jj)) = localptr(ind(jj)) + qf
                                localptrd(ind(jj)) = localptrd(ind(jj)) + qfd
                            end do
                        end do
                    end do
                    ii = ii + ni * nj
                end if
            end do
        end do
 
        call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecrestorearrayf90(nodevallocald, localptrd, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Globalize the current force
        call vecset(exch%nodeValGlobal, zero, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                             exch%nodeValGlobal, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, exch%nodeValLocal, &
                           exch%nodeValGlobal, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Globalize the current force derivative
        call vecset(nodevalglobald, zero, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterbegin(exch%scatter, nodevallocald, &
                             nodevalglobald, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, nodevallocald, &
                           nodevalglobald, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        ! The product rule here: (since we are multiplying)
        ! nodeValGlobal = nodeValGlobal * invArea
        ! nodeValGlobald = nodeValGlobald*invArea + nodeValGlobal*invAread
 
        ! First term:  nodeValGlobald = nodeValGlobald*invArea
        call vecpointwisemult(nodevalglobald, nodevalglobald, &
                              exch%sumGlobal, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Second term:, tmp = nodeValGlobal*invAread
        call vecpointwisemult(tmp, exch%nodeValGlobal, sumglobald, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Sum the second term into the first
        call vecaxpy(nodevalglobald, one, tmp, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Push back to the local values
        call vecscatterbegin(exch%scatter, nodevalglobald, &
                             nodevallocald, insert_values, scatter_reverse, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, nodevalglobald, &
                           nodevallocald, insert_values, scatter_reverse, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecgetarrayf90(nodevallocald, localptrd, ierr)
        call echk(ierr, __file__, __line__)
 
        ii = 0
        do nn = 1, ndom
            call setpointers_d(nn, 1_inttype, sps)
            do mm = 1, nbocos
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
 
                if (bctype(mm) == eulerwall .or. &
                    bctype(mm) == nswalladiabatic .or. &
                    bctype(mm) == nswallisothermal) then
                    do j = jbeg, jend
                        do i = ibeg, iend
                            ii = ii + 1
                            if (idim <= 3) then
                                bcdatad(mm)%Tp(i, j, idim) = localptrd(ii)
                            else
                                bcdatad(mm)%Tv(i, j, idim - 3) = localptrd(ii)
                            end if
                        end do
                    end do
                end if
            end do
        end do
 
        call vecrestorearrayf90(nodevallocald, localptrd, ierr)
        call echk(ierr, __file__, __line__)
 
    end do dimloop
 
    call vecdestroy(nodevallocald, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(nodevalglobald, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(sumglobald, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(tmp, ierr)
    call echk(ierr, __file__, __line__)
 
end subroutine computenodaltractions_d
 
subroutine computenodaltractions_b(sps)
 
    ! This routine performs the reverse of computeNodalTractions. Tt
    ! takes in bcDatad%Tv and bcDatad%Tp and perfroms the reverse of the
    ! nodal averaging procedure in getForces to compute bcDatad(mm)%Fp,
    ! bcDatad(mm)%Fv and bcDatad(mm)%area.
 
    use constants
    use blockpointers, only: ndom, nbocos, bcdata, bctype, nbocos, bcdatad
    use inputphysics, only: forcesastractions
    use surfacefamilies, only: bcfamexchange, familyexchange
    use communication
    use utils, only: echk, setpointers, setpointers_d
#include <petsc/finclude/petsc.h>
    use petsc
    implicit none
 
    integer(kind=intType), intent(in) :: sps
    integer(kind=intType) :: mm, nn, i, j, ii, jj, idim, ierr
    integer(kind=intType) :: ibeg, iend, jbeg, jend, ind(4), ni, nj
    real(kind=realtype) :: qf_b, qf, qa, qa_b
    real(kind=realtype), dimension(:), pointer :: localptr, localptr_b
    type(familyexchange), pointer :: exch
    vec nodevallocal_b, nodevalglobal_b, sumglobal_b, tmp, tmp_b, t_b
 
    ! For better readibility
    exch => bcfamexchange(ibcgroupwalls, sps)
 
    call vecduplicate(exch%nodeValLocal, nodevallocal_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(exch%nodeValGlobal, nodevalglobal_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(exch%sumGlobal, sumglobal_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(exch%sumGlobal, tmp, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecduplicate(tmp, t_b, ierr)
    call echk(ierr, __file__, __line__)
 
    ! For tractions it's (a lot) more difficult becuase we have to do
    ! the scatter/gather operation.
 
    ! ==================================
    !  Recompute the dual area
    ! ==================================
 
    call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    localptr = zero
    ! ii is the running counter through the pointer array.
    ii = 0
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
            ni = iend - ibeg + 1
            nj = jend - jbeg + 1
 
            if (bctype(mm) == eulerwall .or. &
                bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
                do j = 0, nj - 2
                    do i = 0, ni - 2
 
                        ! Scatter a quarter of the area to each node:
                        qa = fourth * bcdata(mm)%area(i + ibeg + 1, j + jbeg + 1)
                        ind(1) = ii + (j) * ni + i + 1
                        ind(2) = ii + (j) * ni + i + 2
                        ind(3) = ii + (j + 1) * ni + i + 2
                        ind(4) = ii + (j + 1) * ni + i + 1
                        do jj = 1, 4
                            localptr(ind(jj)) = localptr(ind(jj)) + qa
                        end do
                    end do
                end do
                ii = ii + ni * nj
            end if
        end do
    end do
 
    call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Globalize the area
    call vecset(exch%sumGlobal, zero, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                         exch%sumGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, exch%nodeValLocal, &
                       exch%sumGlobal, add_values, scatter_forward, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Now compute the inverse of the weighting so that we can multiply
    ! instead of dividing.
 
    call vecgetarrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    localptr = one / localptr
 
    call vecrestorearrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ! ==================================
    ! Now trace through the computeNodalTractions() routine
    ! backwards. All the scatters flip direction and INSERT_VALUES
    ! becomes ADD_VALUES and vice-versa
    ! ==================================
    dimloop: do idim = 1, 6
 
        ! ====================
        ! Do the forward pass:
        ! ====================
        call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        localptr = zero
        ! ii is the running counter through the pointer array.
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, 1_inttype, sps)
            do mm = 1, nbocos
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                if (bctype(mm) == eulerwall .or. &
                    bctype(mm) == nswalladiabatic .or. &
                    bctype(mm) == nswallisothermal) then
                    do j = 0, nj - 2
                        do i = 0, ni - 2
                            if (idim <= 3) then
                                qf = fourth * bcdata(mm)%Fp(i + ibeg + 1, j + jbeg + 1, idim)
                            else
                                qf = fourth * bcdata(mm)%Fv(i + ibeg + 1, j + jbeg + 1, idim - 3)
                            end if
 
                            ind(1) = ii + (j) * ni + i + 1
                            ind(2) = ii + (j) * ni + i + 2
                            ind(3) = ii + (j + 1) * ni + i + 2
                            ind(4) = ii + (j + 1) * ni + i + 1
                            do jj = 1, 4
                                localptr(ind(jj)) = localptr(ind(jj)) + qf
                            end do
                        end do
                    end do
                    ii = ii + ni * nj
                end if
            end do
        end do
 
        call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Globalize the current force
        call vecset(exch%nodeValGlobal, zero, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                             exch%nodeValGlobal, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, exch%nodeValLocal, &
                           exch%nodeValGlobal, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        ! ====================
        ! Do the reverse pass:
        ! ====================
 
        ! Copy the reverse seed into the local values
        call vecgetarrayf90(nodevallocal_b, localptr_b, ierr)
        call echk(ierr, __file__, __line__)
 
        ii = 0
        domains: do nn = 1, ndom
            call setpointers_d(nn, 1_inttype, sps)
 
            ! Loop over the number of boundary subfaces of this block.
            bocos: do mm = 1, nbocos
                if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                    bctype(mm) == nswallisothermal) then
 
                    ! This is easy, just copy out F or T in continuous ordering.
                    do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd
                        do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd
                            ii = ii + 1
                            if (idim <= 3) then
                                localptr_b(ii) = bcdatad(mm)%Tp(i, j, idim)
                            else
                                localptr_b(ii) = bcdatad(mm)%Tv(i, j, idim - 3)
                            end if
                        end do
                    end do
                end if
            end do bocos
        end do domains
 
        call vecrestorearrayf90(nodevallocal_b, localptr_b, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecset(t_b, zero, ierr)
        call echk(ierr, __file__, __line__)
        ! Push up to the global values
        call vecscatterbegin(exch%scatter, nodevallocal_b, &
                             t_b, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, nodevallocal_b, &
                           t_b, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        ! this is particularly nasty. This is why you don't do
        ! derivatives by hand, kids.
        ! exch%nodeValGlobal = F
        ! nodeValGlobal_b = F_b
        ! T_b = reverse seed for tractions
        ! sumGlobal_b =  inverseDualarea_b
        ! exch%sumGlobal = invDualArea
 
        ! Basically what we have to compute here is:
        ! Fb = invDualArea * T_b
        ! invDualAreab = invDualAreab + F*T_b
 
        call vecpointwisemult(nodevalglobal_b, exch%sumGlobal, t_b, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecpointwisemult(tmp, exch%nodeValGlobal, t_b, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Accumulate seed on adflowGlobal_b
        call vecaxpy(sumglobal_b, one, tmp, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Now communicate F_b back to the local patches
 
        call vecscatterbegin(exch%scatter, nodevalglobal_b, &
                             nodevallocal_b, insert_values, scatter_reverse, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, nodevalglobal_b, &
                           nodevallocal_b, insert_values, scatter_reverse, ierr)
        call echk(ierr, __file__, __line__)
 
        ! ============================
        ! Copy the values into patches
        ! ============================
 
        call vecgetarrayf90(nodevallocal_b, localptr_b, ierr)
        call echk(ierr, __file__, __line__)
 
        ! ii is the running counter through the pointer array.
        ii = 0
        do nn = 1, ndom
            call setpointers_d(nn, 1_inttype, sps)
            do mm = 1, nbocos
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                if (bctype(mm) == eulerwall .or. &
                    bctype(mm) == nswalladiabatic .or. &
                    bctype(mm) == nswallisothermal) then
 
                    do j = 0, nj - 2
                        do i = 0, ni - 2
 
                            ind(1) = ii + (j) * ni + i + 1
                            ind(2) = ii + (j) * ni + i + 2
                            ind(3) = ii + (j + 1) * ni + i + 2
                            ind(4) = ii + (j + 1) * ni + i + 1
                            qf_b = zero
                            do jj = 1, 4
                                qf_b = qf_b + localptr_b(ind(jj))
                            end do
                            qf_b = qf_b * fourth
 
                            if (idim <= 3) then
                                bcdatad(mm)%Fp(i + ibeg + 1, j + jbeg + 1, idim) = &
                                    bcdatad(mm)%Fp(i + ibeg + 1, j + jbeg + 1, idim) + qf_b
                            else
                                bcdatad(mm)%Fv(i + ibeg + 1, j + jbeg + 1, idim - 3) = &
                                    bcdatad(mm)%Fv(i + ibeg + 1, j + jbeg + 1, idim - 3) + qf_b
                            end if
 
                        end do
                    end do
                    ii = ii + ni * nj
                end if
            end do
        end do
 
        call vecrestorearrayf90(nodevallocal_b, localptr_b, ierr)
        call echk(ierr, __file__, __line__)
 
    end do dimloop
 
    ! ============================
    ! Finish the dual area sensitivity.
    ! ============================
 
    ! On the forward pass we computed:
    ! sumGlobal = one/sumGlobal
    ! So on the reverse pass we need:
    ! sumGlobalb = -(sumGlobalb/sumGlobal**2)
 
    ! We will do this by getting pointers
 
    call vecgetarrayf90(sumglobal_b, localptr_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecgetarrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Keep in mind localPtr points to sumGlobal which already has
    ! been inversed so we just multiply.
    localptr_b = -localptr_b * localptr**2
 
    call vecrestorearrayf90(sumglobal_b, localptr_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecrestorearrayf90(exch%sumGlobal, localptr, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Push back to the local patches
    call vecscatterbegin(exch%scatter, sumglobal_b, &
                         nodevallocal_b, insert_values, scatter_reverse, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecscatterend(exch%scatter, sumglobal_b, &
                       nodevallocal_b, insert_values, scatter_reverse, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecgetarrayf90(nodevallocal_b, localptr_b, ierr)
    call echk(ierr, __file__, __line__)
 
    ! ii is the running counter through the pointer array.
    ii = 0
    do nn = 1, ndom
        call setpointers_d(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
            ni = iend - ibeg + 1
            nj = jend - jbeg + 1
 
            if (bctype(mm) == eulerwall .or. &
                bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
                do j = 0, nj - 2
                    do i = 0, ni - 2
 
                        ind(1) = ii + (j) * ni + i + 1
                        ind(2) = ii + (j) * ni + i + 2
                        ind(3) = ii + (j + 1) * ni + i + 2
                        ind(4) = ii + (j + 1) * ni + i + 1
                        qa_b = zero
                        do jj = 1, 4
                            qa_b = qa_b + localptr_b(ind(jj))
                        end do
                        qa_b = fourth * qa_b
                        bcdatad(mm)%area(i + ibeg + 1, j + jbeg + 1) = &
                            bcdatad(mm)%area(i + ibeg + 1, j + jbeg + 1) + qa_b
                    end do
                end do
                ii = ii + ni * nj
            end if
        end do
    end do
 
    call vecrestorearrayf90(nodevallocal_b, localptr_b, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Remove temporary petsc vecs
    call vecdestroy(nodevallocal_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(nodevalglobal_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(sumglobal_b, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(tmp, ierr)
    call echk(ierr, __file__, __line__)
 
    call vecdestroy(t_b, ierr)
    call echk(ierr, __file__, __line__)
 
end subroutine computenodaltractions_b
 
subroutine computenodalforces(sps)
 
    ! This subroutine averages the cell based forces and tractions to
    ! node based values. There is no need for communication since we are
    ! simplying summing a quarter of each value to each corner.
 
    use constants
    use blockpointers, only: ndom, nbocos, bctype, bcdata
    use utils, only: setpointers
    implicit none
 
    integer(kind=intType), intent(in) :: sps
 
    integer(kind=intType) :: mm, nn, i, j
    integer(kind=intType) :: ibeg, iend, jbeg, jend
    real(kind=realtype) :: qf(3)
 
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
 
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
                bcdata(mm)%F = zero
                do j = jbeg, jend
                    do i = ibeg, iend
                        qf = fourth * (bcdata(mm)%Fp(i, j, :) + bcdata(mm)%Fv(i, j, :))
                        bcdata(mm)%F(i, j, :) = bcdata(mm)%F(i, j, :) + qf
                        bcdata(mm)%F(i - 1, j, :) = bcdata(mm)%F(i - 1, j, :) + qf
                        bcdata(mm)%F(i, j - 1, :) = bcdata(mm)%F(i, j - 1, :) + qf
                        bcdata(mm)%F(i - 1, j - 1, :) = bcdata(mm)%F(i - 1, j - 1, :) + qf
                    end do
                end do
            end if
        end do
    end do
end subroutine computenodalforces
 
subroutine computenodalforces_d(sps)
 
    ! Forward mode linearization of nodalForces
 
    use constants
    use blockpointers, only: ndom, nbocos, bctype, bcdata, bcdatad
    use utils, only: setpointers
    implicit none
 
    integer(kind=intType), intent(in) :: sps
 
    integer(kind=intType) :: mm, nn, i, j
    integer(kind=intType) :: ibeg, iend, jbeg, jend
    real(kind=realtype) :: qfd(3)
 
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
 
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
                bcdatad(mm)%F = zero
                do j = jbeg, jend
                    do i = ibeg, iend
                        qfd = fourth * (bcdatad(mm)%Fp(i, j, :) + bcdatad(mm)%Fv(i, j, :))
                        bcdatad(mm)%F(i, j, :) = bcdatad(mm)%F(i, j, :) + qfd
                        bcdatad(mm)%F(i - 1, j, :) = bcdatad(mm)%F(i - 1, j, :) + qfd
                        bcdatad(mm)%F(i, j - 1, :) = bcdatad(mm)%F(i, j - 1, :) + qfd
                        bcdatad(mm)%F(i - 1, j - 1, :) = bcdatad(mm)%F(i - 1, j - 1, :) + qfd
                    end do
                end do
            end if
        end do
    end do
end subroutine computenodalforces_d
 
subroutine computenodalforces_b(sps)
 
    ! Reverse mode linearization of nodalForces
 
    use constants
    use blockpointers, only: ndom, nbocos, bctype, bcdata, bcdatad
    use utils, only: setpointers_d
    implicit none
 
    integer(kind=intType), intent(in) :: sps
 
    integer(kind=intType) :: mm, nn, i, j
    integer(kind=intType) :: ibeg, iend, jbeg, jend
    real(kind=realtype) :: qf_b(3)
 
    domains: do nn = 1, ndom
        call setpointers_d(nn, 1_inttype, sps)
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
                do j = jbeg, jend
                    do i = ibeg, iend
                        qf_b = fourth * (bcdatad(mm)%F(i, j, :) + bcdatad(mm)%F(i - 1, j, :) + &
                                         bcdatad(mm)%F(i, j - 1, :) + bcdatad(mm)%F(i - 1, j - 1, :))
 
                        ! Fp and Fv are face-based values
                        bcdatad(mm)%Fp(i, j, :) = bcdatad(mm)%Fp(i, j, :) + qf_b
                        bcdatad(mm)%Fv(i, j, :) = bcdatad(mm)%Fv(i, j, :) + qf_b
                    end do
                end do
                ! this needs to be after the update to be the reverse of the forward mode.
                bcdatad(mm)%F = zero
            end if
        end do
    end do domains
end subroutine computenodalforces_b
 
subroutine getheatflux(hflux, npts, sps)
    use constants
    use blockpointers, only: ndom, nbocos, bctype, bcdata
    use surfacefamilies, only: bcfamexchange, familyexchange, &
                               zerocellval, zeronodeval
    use utils, only: setpointers
    implicit none
    !
    !      Local variables.
    !
    integer(kind=intType), intent(in) :: npts, sps
    real(kind=realtype), intent(out) :: hflux(npts)
 
    integer(kind=intType) :: mm, nn, i, j, ii
    integer(kind=intType) :: ibeg, iend, jbeg, jend
    type(familyexchange), pointer :: exch
 
    exch => bcfamexchange(ibcgroupwalls, sps)
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        call heatfluxes()
 
        do mm = 1, nbocos
            ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
            jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
 
            bocotype1: if (bctype(mm) == nswallisothermal) then
                bcdata(mm)%cellVal => bcdata(mm)%area(:, :)
            else if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic) then
                bcdata(mm)%cellVal => zerocellval
                bcdata(mm)%nodeVal => zeronodeval
            end if bocotype1
        end do
    end do
 
    call computeweighting(exch)
 
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        do mm = 1, nbocos
            bocotype2: if (bctype(mm) == nswallisothermal) then
                bcdata(mm)%cellVal => bcdata(mm)%cellHeatFlux(:, :)
                bcdata(mm)%nodeVal => bcdata(mm)%nodeHeatFlux(:, :)
            end if bocotype2
        end do
    end do
 
    call surfacecellcentertonode(exch)
 
    ! Now extract into the flat array:
    ii = 0
    do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
 
        ! Loop over the number of viscous boundary subfaces of this block.
        ! According to preprocessing/viscSubfaceInfo, visc bocos are numbered
        ! before other bocos. Therefore, mm_nViscBocos == mm_nBocos
        do mm = 1, nbocos
            bocotype3: if (bctype(mm) == nswallisothermal) then
                do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd
                    do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd
                        ii = ii + 1
                        hflux(ii) = bcdata(mm)%nodeHeatFlux(i, j)
                    end do
                end do
                ! Simply put in zeros for the other wall BCs
            else if (bctype(mm) == nswalladiabatic .or. bctype(mm) == eulerwall) then
                do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd
                    do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd
                        ii = ii + 1
                        hflux(ii) = zero
                    end do
                end do
            end if bocotype3
        end do
    end do
end subroutine getheatflux
 
subroutine heatfluxes
    use constants
    use blockpointers, only: bcdata, ndom, nbocos, bctype, bcfaceid, viscsubface
    use bcpointers, only: ssi
    use flowvarrefstate, only: pref, rhoref
    use utils, only: setpointers, setbcpointers
    implicit none
    !
    !      Local variables.
    !
    integer(kind=intType) :: i, j, ii, mm
    real(kind=realtype) :: fact, scaledim
    real(kind=realtype) :: qw, qa
    logical :: heatedSubface
 
    ! Set the actual scaling factor such that ACTUAL heat flux is computed
    ! The factor is determined from stanton number
    scaledim = pref * sqrt(pref / rhoref)
 
    ! Loop over the boundary subfaces of this block.
    bocos: do mm = 1, nbocos
 
        ! Only do this on isoThermalWalls
        if (bctype(mm) == nswallisothermal) then
 
            ! Set a bunch of pointers depending on the face id to make
            ! a generic treatment possible. The routine setBcPointers
            ! is not used, because quite a few other ones are needed.
            call setbcpointers(mm, .true.)
 
            select case (bcfaceid(mm))
            case (imin, jmin, kmin)
                fact = -one
            case (imax, jmax, kmax)
                fact = one
            end select
 
            ! Loop over the quadrilateral faces of the subface. Note that
            ! the nodal range of BCData must be used and not the cell
            ! range, because the latter may include the halo's in i and
            ! j-direction. The offset +1 is there, because inBeg and jnBeg
            ! refer to nodal ranges and not to cell ranges.
            !
            do j = (bcdata(mm)%jnBeg + 1), bcdata(mm)%jnEnd
                do i = (bcdata(mm)%inBeg + 1), bcdata(mm)%inEnd
 
                    ! Compute the normal heat flux on the face. Inward positive.
                    bcdata(mm)%cellHeatFlux(i, j) = -fact * scaledim * &
                                                    sqrt(ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2) * &
                                                    (viscsubface(mm)%q(i, j, 1) * bcdata(mm)%norm(i, j, 1) &
                                                     + viscsubface(mm)%q(i, j, 2) * bcdata(mm)%norm(i, j, 2) &
                                                     + viscsubface(mm)%q(i, j, 3) * bcdata(mm)%norm(i, j, 3))
                end do
            end do
        end if
    end do bocos
end subroutine heatfluxes
 
subroutine settnswall(tnsw, npts, sps)
 
    use constants
    use blockpointers, only: ndom, nbocos, bcdata, bctype
    use flowvarrefstate, only: tref
    use utils, only: setpointers
    implicit none
 
    ! Input Variables
    integer(kind=intType), intent(in) :: npts, sps
    real(kind=realtype), intent(in) :: tnsw(npts)
 
    ! Local Variables
    integer(kind=intType) :: mm, nn, i, j, ii
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
 
    ii = 0
    domains: do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        ! Loop over the number of viscous boundary subfaces of this block.
        bocos: do mm = 1, nbocos
            isowall: if (bctype(mm) == nswallisothermal) then
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                do j = jbeg, jend
                    do i = ibeg, iend
                        ii = ii + 1
                        bcdata(mm)%TNS_Wall(i, j) = tnsw(ii) / tref
                    end do
                end do
            end if isowall
        end do bocos
    end do domains
    ! TODO: The temperature must be interpolated to the coarse meshes.
    !
    ! The following lines are extracted from BCData/setBCDataCoarseGrid
    ! By the design of the subroutine, TNSWall shall be interpolated during this process.
    ! Yet, the subroutine requires an internal subroutine interpolateBcData.
    ! It remains a question whether interpolateBcData shall become a normal subroutine.
    !
    ! use blockPointers, only : flowDoms
    ! use inputTimeSpectral, only : nTimeIntervalsSpectral
    ! use iteration, only : groundLevel
    ! implicit none
    ! !
    ! !      Local variables.
    ! !
    ! integer(kind=intType) :: nLevels, level, levm1
 
    ! ! Determine the number of grid levels.
 
    ! nLevels = ubound(flowDoms,2)
 
    ! ! Loop over the coarser grid levels. It is assumed that the
    ! ! bc data of groundLevel is set correctly.
 
    ! coarseLevelLoop: do level=(groundLevel+1),nLevels
 
    !    ! Store the fine grid level a bit easier.
 
    !    levm1 = level - 1
 
    !    ! Loop over the number of spectral solutions and local blocks.
 
    !    spectralLoop: do sps=1,nTimeIntervalsSpectral
    !       domainsLoop: do i=1,nDom
 
    !          ! Set the pointers to the coarse block.
 
    !          call setPointers(i, level, sps)
 
    !          ! Loop over the boundary subfaces and interpolate the
    !          ! prescribed boundary data for this grid level.
 
    !          bocoLoop: do j=1,nBocos
 
    !             ! Interpolate the data for the possible prescribed boundary
    !             ! data.
 
    !             call interpolateBcData(BCData(j)%TNS_Wall, &
    !                  flowDoms(i,levm1,sps)%BCData(j)%TNS_Wall)
 
    !          enddo bocoLoop
    !       enddo domainsLoop
    !    enddo spectralLoop
    ! enddo coarseLevelLoop
 
end subroutine settnswall
 
subroutine gettnswall(tnsw, npts, sps)
 
    use constants
    use blockpointers, only: ndom, nbocos, bcdata, bctype
    use flowvarrefstate, only: tref
    use utils, only: setpointers
    implicit none
 
    ! Input Variables
    integer(kind=intType), intent(in) :: npts, sps
    real(kind=realtype), intent(out) :: tnsw(npts)
 
    ! Local Variables
    integer(kind=intType) :: mm, nn, i, j, ii
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
 
    ii = 0
    domains: do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        ! Loop over the number of viscous boundary subfaces of this block.
        bocos: do mm = 1, nbocos
            isowall: if (bctype(mm) == nswallisothermal) then
                do j = jbeg, jend
                    do i = ibeg, iend
                        ii = ii + 1
                        tnsw(ii) = bcdata(mm)%TNS_Wall(i, j) * tref
                    end do
                end do
            end if isowall
        end do bocos
    end do domains
end subroutine gettnswall