zipperMesh.F90

Go to the documentation of this file.

module zippermesh
contains
 
    !
    !      createZipperMesh zips multiple overlapping surface meshes.
    !      First, it eliminates overlapping quads and then stiches the
    !      non-overlapping surface mesh boundaries with triangular
    !      surface grids.
    !      In an overset framework the overlapping surface grids give
    !      wrong estimate of airloads. Zipper mesh overcomes this by
    !      creating a more accurate representation of the surface in the
    !      mesh overlap regions.
    !      ref:  "Enhancements to the Hybrid Mesh Approach to Surface Loads
    !             Integration on Overset Structured Grids", William M. Chan
    !      http://people.nas.nasa.gov/~wchan/publications/AIAA-2009-3990.pdf
    !
 
    subroutine createzippermesh(zipperFamList, nZipFam)
 
        use constants
        use communication, only: myid, adflow_comm_world, nproc, recvrequests, &
                                 sendrequests, commpatterncell_2nd, internalcell_2nd
        use blockpointers, only: ndom, bcdata, nbocos, bctype, il, jl, kl
        use oversetdata, only: oversetstring, oversetwall, csrmatrix, cumdomproc, ndomtotal, &
                               clusters, overlapmatrix, &
                               oversetpresent, zippermesh, zippermeshes
        use walldistancedata, only: xvolumevec, is1, is2
        use utils, only: setpointers, echk
        use adjointvars, only: nnodeslocal
        use sorting, only: faminlist
        use oversetutilities, only: getworkarray, deallocateosurfs, transposeoverlap
        use oversetcommutilities, only: exchangesurfaceiblanks, recvosurf, sendosurf, &
                                        getosurfcommpattern
        use oversetpackingroutines, only: packosurf, unpackosurf, getosurfbuffersizes
        use oversetinitialization, only: initializeosurf
        use inputoverset, only: debugzipper, usezippermesh
        use surfacefamilies, only: bcfamexchange, famnames, bcfamgroups
        use stringops
        use gapboundaries
        use wallsearches, only: wallsearch
 
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: nZipFam
        integer(kind=intType), intent(in), dimension(nZipFam) :: zipperFamList
 
        ! Local Variables
        integer(kind=intType) :: i, j, k, ii, jj, kk, iStart, iSize, sps, level, iStr
        integer(kind=intType) :: iDom, jDom, nn, mm, n, ierr, iProc, iWork, nodeFam(3)
        integer(kind=intType) :: nNodeTotal, nTriTotal, offset, iBCGroup, nFam
        integer(kind=intType), dimension(:), allocatable :: famList
 
        integer(kind=intType), dimension(:, :), allocatable :: tmpInt2D, work
        logical, dimension(:), allocatable :: oSurfReady
        type(oversetwall), dimension(:), allocatable :: oSurfs
        integer(kind=intType), dimension(:), allocatable :: intRecvBuf
        type(oversetstring), target :: master, pocketMaster
        type(oversetstring), pointer :: str
        integer(kind=intType), dimension(:), allocatable :: nodeIndices
 
        integer(kind=intType) :: nFullStrings, nUnique
        integer(kind=intType) :: nOSurfRecv, nOSurfSend
        integer(kind=intType), dimension(:, :), allocatable :: oSurfRecvList, oSurfSendList
        ! MPI/Communication related
        integer mpiStatus(MPI_STATUS_SIZE)
        integer(kind=intType), dimension(:, :), allocatable :: bufSizes
        integer(kind=intType), dimension(:, :), allocatable :: recvInfo
        integer(kind=intType) :: sendCount, recvCount, index
        type(csrmatrix), pointer :: overlap
        type(csrmatrix) :: overlapTranspose
        type(zippermesh), pointer :: zipper
        ! Wall search related
        integer(kind=intType) :: ncells
        type(oversetwall), dimension(:), allocatable, target :: walls
        type(oversetwall), target :: fullWall
 
        if (.not. oversetpresent .or. (.not. usezippermesh)) then
            ! Not overset so we don't can't have a zipper.
            return
        end if
 
        ! Zipper is only implemented for single grid and 1 spectral
        ! instance (ie not time spectral).
        level = 1
        sps = 1
 
        call initbcdataiblank(level, sps)
 
        ! We build the zipper meshes *independently* for each BCType.
        bcgrouploop: do ibcgroup = 1, nfamexchange
 
            ! Set a pointer to the zipper we are working on to make code
            ! easier to read
            zipper => zippermeshes(ibcgroup)
 
            ! Deallocate if already exists
            if (zipper%allocated) then
                call vecscatterdestroy(zipper%scatter, ierr)
                call echk(ierr, __file__, __line__)
                call vecdestroy(zipper%localVal, ierr)
                call echk(ierr, __file__, __line__)
                zipper%allocated = .false.
            end if
 
            ! Note that the zipper%conn could be allocated with size of
            ! zero, but the vecScatter and Vec are not petsc-allocated.
            if (allocated(zipper%conn)) then
                deallocate (zipper%conn, zipper%fam, zipper%indices)
            end if
 
            ! Before we can proceed with the zipper, we need to generate
            ! intersection of the zipperFamList() with the families on this
            ! BCGroup. This would be so much easier in Python...
 
            if (allocated(famlist)) then
                deallocate (famlist)
            end if
 
            nfam = 0
            do i = 1, size(bcfamgroups(ibcgroup)%famList)
                do j = 1, size(zipperfamlist)
                    if (bcfamgroups(ibcgroup)%famList(i) == zipperfamlist(j)) then
                        nfam = nfam + 1
                    end if
                end do
            end do
 
            ! If nFam is zero, no need ot do anything for this zipper. Just
            ! allocated zero-sized arrays so we know the size is 0.
            if (nfam == 0) then
                allocate (zipper%conn(3, 0), zipper%fam(0), zipper%indices(0))
 
                cycle ! to the next BCGroup
            else
                ! Do a second pass and fill up the famList
                allocate (famlist(nfam))
                nfam = 0
                do i = 1, size(bcfamgroups(ibcgroup)%famList)
                    do j = 1, size(zipperfamlist)
                        if (bcfamgroups(ibcgroup)%famList(i) == zipperfamlist(j)) then
                            nfam = nfam + 1
                            famlist(nfam) = zipperfamlist(j)
                        end if
                    end do
                end do
            end if
 
            if (debugzipper .and. myid == 0) then
                write (*, "(a)", advance="no") '-> Creating zipper for families : '
                do i = 1, size(famlist)
                    write (*, "(a,1x)", advance="no") trim(famnames(famlist(i)))
                end do
                print "(1x)"
            end if
 
            overlap => overlapmatrix(level, sps)
            call transposeoverlap(overlap, overlaptranspose)
            ! -------------------------------------------------------------------
            ! Step 1: Eliminate any gap overlaps between meshes
            ! -------------------------------------------------------------------
 
            ! Determine the average area of surfaces on each cluster. This
            ! will be independent of any block splitting distribution.
            call determineclusterareas(famlist)
 
            ! Set the boundary condition blank values
            call slitelimination(famlist, level, sps)
 
            ! Create the surface deltas
            call surfacedeviation(famlist, level, sps)
 
            ! ================ WE NORMALLY GOT THIS FROM OVERSETCOMM =============
            ! Get the OSurf buffer sizes becuase we need that for getOSurfCommPattern
            allocate (bufsizes(ndomtotal, 6), tmpint2d(ndomtotal, 6))
            tmpint2d = 0
            do nn = 1, ndom
                call setpointers(nn, level, sps)
                idom = cumdomproc(myid) + nn
                call getosurfbuffersizes(famlist, il, jl, kl, tmpint2d(idom, 5), &
                                         tmpint2d(idom, 6), .true.)
            end do
 
            call mpi_allreduce(tmpint2d, bufsizes, 6 * ndomtotal, adflow_integer, mpi_sum, &
                               adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Get the basic surface comm pattern.
            ! For sending, the worse case is sending all my blocks/fringes/walls to
            ! everyone but myself:
            ii = ndom * (nproc - 1)
            allocate (osurfsendlist(2, ii))
 
            ! For receiving, the worse receive is all the blocks/fringes/wall I
            ! don't already have:
            ii = ndomtotal - ndom
            allocate (osurfrecvlist(2, ii))
 
            call getosurfcommpattern(overlap, overlaptranspose, &
                                     osurfsendlist, nosurfsend, osurfrecvlist, nosurfrecv, bufsizes(:, 6))
            deallocate (tmpint2d, bufsizes)
            ! ========================================================================
 
            ! Alloc data for the OSurfs
            nn = max(nproc, 2 * nosurfsend, 2 * nosurfrecv)
            allocate (tmpint2d(ndomtotal, 2), bufsizes(ndomtotal, 2), &
                      osurfready(ndomtotal), recvinfo(2, nn))
 
            tmpint2d = 0
            ! Need to get the differt sizes for the OSurfs since they are now
            ! based on the primal mesh as opposed to do the dual mesh as
            ! previously done
            do nn = 1, ndom
                call setpointers(nn, level, sps)
                idom = cumdomproc(myid) + nn
                call getosurfbuffersizes(famlist, il, jl, kl, tmpint2d(idom, 1), &
                                         tmpint2d(idom, 2), .false.)
            end do
 
            ! Make sure everyone has the sizes
            call mpi_allreduce(tmpint2d, bufsizes, 2 * ndomtotal, adflow_integer, mpi_sum, &
                               adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
            deallocate (tmpint2d)
 
            ! allocate oSurfs for the primal mesh
            allocate (osurfs(ndomtotal))
            do idom = 1, ndomtotal
                if (bufsizes(idom, 1) == 0) then
                    osurfready(idom) = .true.
                else
                    osurfready(idom) = .false.
                end if
            end do
 
            ! Initialize the primal walls
            do nn = 1, ndom
                call setpointers(nn, level, sps)
                idom = cumdomproc(myid) + nn
                call initializeosurf(famlist, osurfs(idom), .false., clusters(idom))
                call packosurf(famlist, osurfs(idom), .false.)
                osurfready(idom) = .true.
            end do
 
            ! Post all the oSurf iSends
            sendcount = 0
            do jj = 1, nosurfsend
                iproc = osurfsendlist(1, jj)
                idom = osurfsendlist(2, jj)
                call sendosurf(osurfs(idom), idom, iproc, 0, sendcount)
            end do
 
            recvcount = 0
            do jj = 1, nosurfrecv
                iproc = osurfrecvlist(1, jj)
                idom = osurfrecvlist(2, jj)
                call recvosurf(osurfs(idom), idom, iproc, 0, &
                               bufsizes(idom, 1), bufsizes(idom, 2), recvcount, recvinfo)
            end do
 
            ! Complete all the recives and sends
            do i = 1, recvcount
                call mpi_waitany(recvcount, recvrequests, index, mpistatus, ierr)
                call echk(ierr, __file__, __line__)
            end do
 
            do i = 1, sendcount
                call mpi_waitany(sendcount, sendrequests, index, mpistatus, ierr)
                call echk(ierr, __file__, __line__)
            end do
 
            ! Unpack any blocks we received if necessary:
            do i = 1, recvcount
 
                ! Global domain index of the recv that finished
                idom = recvinfo(1, i)
                if (.not. osurfs(idom)%allocated) then
                    call unpackosurf(osurfs(idom))
                end if
            end do
 
            ! Determine the size of the buffer we need locally for the
            ! receives. We do this outside the main iteration loop since it
            ! always the same size.
            ii = 0
            do jj = 1, nosurfsend
                ! These blocks are by definition local.
                idom = osurfsendlist(2, jj)
                ii = ii + osurfs(idom)%maxCells
            end do
            allocate (intrecvbuf(max(1, ii)))
 
            ! ------------------------ Performing Searches ----------------
            call getworkarray(overlap, work)
 
            do iwork = 1, size(work, 2)
                idom = work(1, iwork)
                jdom = work(2, iwork)
                call wallsearch(osurfs(idom), osurfs(jdom))
            end do
 
            ! ------------------------ Receiving iBlank back ---------------
 
            sendcount = 0
            do jj = 1, nosurfrecv
                iproc = osurfrecvlist(1, jj)
                idom = osurfrecvlist(2, jj)
                sendcount = sendcount + 1
                call mpi_isend(osurfs(idom)%iBlank, osurfs(idom)%maxCells, adflow_integer, &
                               iproc, idom, adflow_comm_world, sendrequests(sendcount), ierr)
                call echk(ierr, __file__, __line__)
            end do
 
            recvcount = 0
            istart = 1
            do jj = 1, nosurfsend
                iproc = osurfsendlist(1, jj)
                idom = osurfsendlist(2, jj)
                isize = osurfs(idom)%maxCells
                recvcount = recvcount + 1
 
                call mpi_irecv(intrecvbuf(istart), isize, adflow_integer, &
                               iproc, idom, adflow_comm_world, recvrequests(recvcount), ierr)
                call echk(ierr, __file__, __line__)
                istart = istart + isize
            end do
 
            ! Now wait for the sends and receives to finish
            do i = 1, sendcount
                call mpi_waitany(sendcount, sendrequests, index, mpistatus, ierr)
                call echk(ierr, __file__, __line__)
            end do
 
            do i = 1, recvcount
                call mpi_waitany(recvcount, recvrequests, index, mpistatus, ierr)
                call echk(ierr, __file__, __line__)
            end do
 
            ! Process the oSurfs we own locally
            do nn = 1, ndom
                call setpointers(nn, level, sps)
                idom = cumdomproc(myid) + nn
                ii = 0
                do mm = 1, nbocos
                    if (faminlist(bcdata(mm)%famID, famlist)) then
                        do j = bcdata(mm)%jnBeg + 1, bcdata(mm)%jnEnd
                            do i = bcdata(mm)%inBeg + 1, bcdata(mm)%inEnd
                                ii = ii + 1
                                if (osurfs(idom)%iBlank(ii) == -2) then
                                    bcdata(mm)%iblank(i, j) = -2
                                end if
                            end do
                        end do
                    end if
                end do
            end do
 
            ! And update based on the data we received from other processors
            ii = 0
            do kk = 1, nosurfsend
 
                idom = osurfsendlist(2, kk)
                nn = idom - cumdomproc(myid)
 
                ! Set the block pointers for the local block we are dealing
                !with:
                call setpointers(nn, level, sps)
                do mm = 1, nbocos
                    if (faminlist(bcdata(mm)%famID, famlist)) then
                        do j = bcdata(mm)%jnBeg + 1, bcdata(mm)%jnEnd
                            do i = bcdata(mm)%inBeg + 1, bcdata(mm)%inEnd
                                ii = ii + 1
                                if (intrecvbuf(ii) == -2) then
                                    bcdata(mm)%iBlank(i, j) = -2
                                end if
                            end do
                        end do
                    end if
                end do
            end do
            ! Release some unnecessary memory
            deallocate (bufsizes, osurfsendlist, osurfrecvlist, osurfready, recvinfo, work)
 
            ! Ditch our oSurfs
            call deallocateosurfs(osurfs, ndomtotal)
            deallocate (osurfs, intrecvbuf)
 
            ! Before we continue, we do a little more
            ! processing. bowTieElimination tries to eliminate cells
            call bowtieandisolationelimination(famlist, level, sps)
 
            ! -------------------------------------------------------------------
            ! Step 2: Identify gap boundary strings and split the strings to
            !         sub-strings.
            ! -------------------------------------------------------------------
 
            if (debugzipper) then
                call writewalls(famlist)
            end if
 
            call makegapboundarystrings(famlist, level, sps, master)
 
            rootproc: if (myid == 0) then
                if (debugzipper) then
                    call writezipperdebug(master)
                end if
 
                ! Run the common core zipper routines
                call zippercore(master, pocketmaster, debugzipper)
 
                ! Before we create the final data structures for the zipper
                ! mesh; we will combine the master and pocketMaster
                ! strings, and unique-ify the indices to help keep amount of data
                ! transfer to a minimum. We also determine at this point which
                ! triangle belongs to which family.
 
                call setstringpointers(master)
                ntritotal = master%ntris + pocketmaster%ntris
                nnodetotal = master%nNodes + pocketmaster%nNodes
                allocate (zipper%conn(3, ntritotal), zipper%fam(ntritotal), zipper%indices(nnodetotal))
 
                ii = 0
                jj = 0
                outerloop: do istr = 1, 2
                    ! Select which of the two we are dealing with
                    if (istr == 1) then
                        str => master
                        offset = 0
                    else
                        str => pocketmaster
                        offset = master%nNodes
                    end if
 
                    ! Loop over the number of triangles.
                    do i = 1, str%nTris
 
                        ! Extract the family from the nodes
                        do j = 1, 3
                            nodefam(j) = str%family(str%tris(j, i))
                        end do
 
                        ! Increment the running counter for all triangles.
                        ii = ii + 1
 
                        ! Set the family
                        zipper%Fam(ii) = selectnodefamily(nodefam)
 
                        ! Set the connectivity.
                        zipper%conn(:, ii) = str%tris(:, i) + offset
                    end do
 
                    ! Loop over the nodal indices and add
                    do j = 1, str%nNodes
                        ! And set the global indices that the zipper needs. Note
                        ! that we are doing zipperIndices as a scalar and that
                        ! the indices refer to the global nodes so are already in
                        ! 0-based ordering.
                        jj = jj + 1
                        zipper%indices(jj) = str%ind(j)
                    end do
                end do outerloop
 
                call deallocatestring(master)
                call deallocatestring(pocketmaster)
            else
                ! Other procs don't have any triangles *sniffle* :-(
                allocate (zipper%conn(3, 0), zipper%fam(0), zipper%indices(0))
            end if rootproc
 
            call veccreatempi(adflow_comm_world, size(zipper%indices), petsc_determine, &
                              zipper%localVal, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Now create the general scatter that goes from the
            ! globalNodalVector to the local vectors.
            call iscreategeneral(adflow_comm_world, size(zipper%indices), &
                                 zipper%indices - 1, petsc_copy_values, is1, ierr)
            call echk(ierr, __file__, __line__)
 
#if PETSC_VERSION_GE(3,8,0)
            call vecscattercreate(bcfamexchange(ibcgroup, sps)%nodeValLocal, is1, &
                                  zipper%localVal, petsc_null_is, zipper%scatter, ierr)
#else
            call vecscattercreate(bcfamexchange(ibcgroup, sps)%nodeValLocal, is1, &
                                  zipper%localVal, petsc_null_object, zipper%scatter, ierr)
#endif
            call echk(ierr, __file__, __line__)
 
            call isdestroy(is1, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Flag to keep track of allocated PETSc object.
            zipper%allocated = .true.
 
        end do bcgrouploop
 
    contains
        function selectnodefamily(nodeFam)
            implicit none
            integer(kind=intType), dimension(3) :: nodefam
            integer(kind=intType) :: selectnodefamily
 
            ! We have a few cases to check:
 
            if (nodefam(1) == nodefam(2) .and. nodefam(1) == nodefam(3)) then
                ! Case 1: All of the nodes have the same family. This is what
                ! *should* happen all the time:
 
                selectnodefamily = nodefam(1)
 
            else if (nodefam(1) == nodefam(2)) then
 
                ! Case 2a: First two nodes are the same. Take the family from 1.
 
                selectnodefamily = nodefam(1)
 
            else if (nodefam(2) == nodefam(3)) then
 
                ! Case 2b: Last two nodes are the same. Take the family from 2.
 
                selectnodefamily = nodefam(2)
 
            else if (nodefam(1) == nodefam(3)) then
 
                ! Case 2b: First and last nodes are the same. Take the family from 1.
 
                selectnodefamily = nodefam(1)
 
            else
 
                ! All nodes are different. We arbitrarily take the first and
                ! print a warning becuase this should not happen.
 
                selectnodefamily = nodefam(1)
                print *, 'Family for triangle could not be uniquely determined. Nodes are from 3 different families!'
            end if
 
        end function selectnodefamily
    end subroutine createzippermesh
 
    subroutine checkzipper(fileName)
 
        ! Special routine for checking zipper mesh loaded from debug file
        use constants
        use stringops
        implicit none
 
        ! Input/Output
        character(*), intent(in) :: fileName
 
        ! Working
        type(oversetstring) :: master, pocketMaster
 
        call loadzipperdebug(filename, master)
        call zippercore(master, pocketmaster, .true.)
 
        print *, 'Zipper successfully completed'
    end subroutine checkzipper
 
    subroutine zippercore(master, pocketMaster, debugZipper)
 
        ! Common routine for creating the zipper from a given master
        use constants
        use stringops
        use kdtree2_module
        implicit none
 
        ! Input/Output
        type(oversetstring), intent(inout) :: master, pocketMaster
        logical, intent(in) :: debugZipper
 
        ! Local
        type(oversetstring), dimension(:), allocatable, target :: strings
        type(oversetstring), pointer :: str
        integer(kind=intType) :: nStrings, i, j, nTriSelf
 
        if (debugzipper) then
            open (unit=101, file="master_beforeStrings.dat", form='formatted')
            write (101, *) 'TITLE = "Master Data" '
            write (101, *) 'Variables = "X" "Y" "Z"'
            call writeoversetmaster(master, 101)
            close (101)
        end if
 
        call createorderedstrings(master, strings, nstrings)
 
        master%myID = 99
 
        ! Allocate space for the maximum number of directed edges. This
        ! is equal to the initial number of edges (nElems) plus 3 times
        ! the number of triangles we will add, which is also nElems. Now,
        ! we will probably not actualy have that many since we save a
        ! triangle and 3 edges for every self zip that is
        ! applied. Therefore we know this will always be enough
        allocate (master%edges(4 * master%nElems))
 
        master%nEdges = 0
 
        do i = 1, nstrings
            str => strings(i)
            do j = 1, str%nElems
                master%nEdges = master%nEdges + 1
                master%edges(master%nEdges)%n1 = str%p%conn(1, str%pElems(j)) 
                master%edges(master%nEdges)%n2 = str%p%conn(2, str%pElems(j)) 
            end do
        end do
 
        ! Allocate space for the triangles. Again, this can be at most,
        ! nElems, but the total number of elements will most likely be
        ! smaller due to self zipping. If someone puts
        allocate (master%tris(3, master%nElems))
        master%nTris = 0
 
        ! Build the master tree
        master%tree => kdtree2_create(master%x, sort=.true.)
 
        ! Perform the string association:
        call stringmatch(strings, nstrings, debugzipper)
 
        if (debugzipper) then
            open (unit=101, file="strings_beforeSelfZip.dat", form='formatted')
            write (101, *) 'TITLE = "Gap Strings Data" '
            write (101, *) 'Variables = "X" "Y" "Z" "Nx" "Ny" "Nz" "Vx" "Vy" "Vz" "ind" &
                 &"gapID" "gapIndex" "otherID" "otherIndex" "ratio"'
            do i = 1, nstrings
                call writeoversetstring(strings(i), strings, nstrings, 101)
            end do
            close (101)
        end if
 
        call performselfzip(master, strings, nstrings, debugzipper)
 
        ! Write out any self-zipped triangles
        ntriself = master%nTris
 
        if (debugzipper) then
            call writeoversettriangles(master, "selfzipTriangulation.dat", 1, master%nTris)
        end if
 
        ! Write out the gaps AFTER the self zip
        if (debugzipper) then
            open (unit=101, file="strings_afterSelfZip.dat", form='formatted')
            write (101, *) 'TITLE = "Gap Strings Data" '
            write (101, *) 'Variables = "X" "Y" "Z" "Nx" "Ny" "Nz" "Vx" "Vy" "Vz" "ind" &
                 &"gapID" "gapIndex" "otherID" "otherIndex" "ratio"'
            do i = 1, nstrings
                call writeoversetstring(strings(i), strings, nstrings, 101)
            end do
            close (101)
        end if
 
        ! Now do the cross zip
        call makecrosszip(master, strings, nstrings, debugzipper)
 
        ! And write out the triangle from the cross zip
        if (debugzipper) then
            call writeoversettriangles(master, "crossZipTriangulation.dat", ntriself + 1, master%nTris)
        end if
 
        ! ---------------------------------------------------------------
        ! Sort through zipped triangle edges and the edges which have not
        ! been used twice (orphan edges) will be ultimately gathered to
        ! form polygon pockets to be zipped.
        if (debugzipper) then
            print *, 'Doing pocket zip'
        end if
        call makepocketzip(master, strings, nstrings, pocketmaster, debugzipper)
 
        if (debugzipper) then
            call writeoversettriangles(pocketmaster, "pocketTriangulation.dat", 1, pocketmaster%nTris)
        end if
 
        ! Clean up the reminder of the sting memory on the root proc
        do i = 1, nstrings
            call deallocatestring(strings(i))
        end do
        deallocate (strings)
 
    end subroutine zippercore
 
    !
    !       determineClusterArea determine the average cell surface area
    !       for all blocks in a particular cluster. This is used for
    !       determine blanking preference for overlapping surface cells.
 
    subroutine determineclusterareas(famList)
 
        use constants
        use blockpointers, only: ndom, bcdata, nbocos, bctype
        use communication, only: adflow_comm_world, myid
        use oversetdata, onlY: clusterareas, nclusters, clusters, cumdomproc
        use utils, only: setpointers, echk, setbcpointers, cross_prod
        use bcpointers, only: xx
        use sorting, only: faminlist
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in), dimension(:) :: famList
 
        ! Working
        integer(kind=intType) :: i, j, mm, nn, clusterID, ierr, nPts, nCells
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
        real(kind=realtype), dimension(:), allocatable :: localareas
        integer(kind=intType), dimension(:), allocatable :: localCount, globalCount
        real(kind=realtype), dimension(:, :), allocatable :: pts
        real(kind=realtype) :: fact, v1(3), v2(3), sss(3), da
 
        if (allocated(clusterareas)) then
            ! We only ever do this once!
            deallocate (clusterareas)
        end if
 
        allocate (clusterareas(nclusters), localareas(nclusters), &
                  localcount(nclusters), globalcount(nclusters))
 
        localareas = zero
        localcount = 0
 
        domains: do nn = 1, ndom
            call setpointers(nn, 1_inttype, 1_inttype)
 
            clusterid = clusters(cumdomproc(myid) + nn)
 
            ! Loop over the number of boundary subfaces of this block.
            bocos: do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
                    ! Store the cell range of the subfaces a bit easier.
                    ! As only owned faces must be considered the nodal range
                    ! in BCData must be used to obtain this data.
 
                    jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                    call setbcpointers(mm, .true.)
 
                    ! Compute the dual area at each node. Just store in first dof
                    do j = jbeg, jend ! This is a face loop
                        do i = ibeg, iend ! This is a face loop
 
                            v1(:) = xx(i + 1, j + 1, :) - xx(i, j, :)
                            v2(:) = xx(i, j + 1, :) - xx(i + 1, j, :)
 
                            ! Cross Product
                            call cross_prod(v1, v2, sss)
                            da = half * sqrt(sss(1)**2 + sss(2)**2 + sss(3)**2)
                            localareas(clusterid) = localareas(clusterid) + da
                            localcount(clusterid) = localcount(clusterid) + 1
                        end do
                    end do
                end if faminclude
            end do bocos
        end do domains
 
        ! All reduce sum for the localAreas to get clusterAreas and
        ! localCount to get globalCount
 
        call mpi_allreduce(localareas, clusterareas, nclusters, adflow_real, &
                           mpi_sum, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        call mpi_allreduce(localcount, globalcount, nclusters, adflow_integer, &
                           mpi_sum, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Final get the average global area
        do i = 1, nclusters
            clusterareas(i) = clusterareas(i) / max(globalcount(i), 1)
        end do
 
        deallocate (localareas, localcount, globalcount)
    end subroutine determineclusterareas
 
    subroutine initbcdataiblank(level, sps)
 
        use constants
        use blockpointers
        use communication
        use utils, only: setpointers
        use oversetutilities, only: flagforcedrecv
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: level, sps
 
        ! Local variables
        integer(kind=intType) :: mm, nn, i, j, k, iBeg, iEnd, jBeg, jEnd
        logical :: side(4)
 
        integer(kind=intType), dimension(:, :), pointer :: ibp, gcp, frx
        integer(kind=intType), dimension(:, :), allocatable :: toFlip
 
        ! This routine initializes the surface cell iblank based on the
        ! volume iblank.
 
        domainloop: do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Setting the surface IBlank array is done for *all* bocos.
            bocoloop: do mm = 1, nbocos
                select case (bcfaceid(mm))
                case (imin)
                    ibp => iblank(2, :, :)
                case (imax)
                    ibp => iblank(il, :, :)
                case (jmin)
                    ibp => iblank(:, 2, :)
                case (jmax)
                    ibp => iblank(:, jl, :)
                case (kmin)
                    ibp => iblank(:, :, 2)
                case (kmax)
                    ibp => iblank(:, :, kl)
                end select
 
                ! -------------------------------------------------
                ! Step 1: Set the (haloed) cell iBlanks directly from
                ! the volume iBlanks
                ! -------------------------------------------------
                jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                ! Just set the cell iblank directly from the cell iblank
                ! above it. Remember the +1 in ibp is for the pointer
                ! offset. These ranges *ALWAYS* give 1 level of halos
                do j = jbeg - 1, jend + 1
                    do i = ibeg - 1, iend + 1
                        bcdata(mm)%iBlank(i, j) = ibp(i + 1, j + 1)
                    end do
                end do
            end do bocoloop
        end do domainloop
    end subroutine initbcdataiblank
 
    subroutine slitelimination(famList, level, sps)
 
        use constants
        use blockpointers
        use communication
        use utils, only: setpointers
        use oversetutilities, only: flagforcedrecv
        use sorting, only: faminlist
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in), dimension(:) :: famList
        integer(kind=intType), intent(in) :: level, sps
 
        ! Local variables
        integer(kind=intType) :: mm, nn, i, j, k, iBeg, iEnd, jBeg, jEnd
        logical :: side(4)
 
        integer(kind=intType), dimension(:, :), pointer :: ibp, gcp, frx
        integer(kind=intType), dimension(:, :), allocatable :: toFlip
 
        ! This routine initializes the surface cell iblank based on the
        ! volume iblank. It is not a straight copy since we a little
        ! preprocessing
 
        ! This is a little trickier than it seems. The reason is that we
        ! will allow a limited number of interpolated cells to be used
        ! directly in the integration provided the meet certain criteria.
 
        ! Consider the following
        !  +------+--------+-------+
        !  |ib=1  |  ib=1  | ib= 1 |
        !  |      |        |       |
        !  |      |        |       |
        !  +------+========+-------+
        !  |ib=1  || ib=-1 || ib=1 |
        !  |      ||       ||      |
        !  |      ||       ||      |
        !==+======+--------+=======+==
        !  |ib=-1 |  ib=-1 | ib=-1 |
        !  |      |        |       |
        !  |      |        |       |
        !  +------+--------+-------+
        !
        ! The boundary between real/interpolated cells is marked by double
        ! lines. For zipper mesh purposes, it is generally going to be
        ! better to treat the center cell, as a regular force integration
        ! cell (ie surface iblank=1). The criteria for selection of these
        ! cells is:
 
        ! 1. The cell must not have been a forced receiver (ie at overset outer
        ! bound)
        ! 2. Any pair *opposite* sides of the cell must be compute cells.
 
        ! This criterial allows one-cell wide 'slits' to be pre-eliminated.
 
        call flagforcedrecv()
        domainloop: do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Setting the surface IBlank array is done for *all* bocos.
            bocoloop: do mm = 1, nbocos
                select case (bcfaceid(mm))
                case (imin)
                    ibp => iblank(2, :, :)
                    gcp => globalcell(2, :, :)
                    frx => forcedrecv(2, :, :)
                case (imax)
                    ibp => iblank(il, :, :)
                    gcp => globalcell(il, :, :)
                    frx => forcedrecv(il, :, :)
                case (jmin)
                    ibp => iblank(:, 2, :)
                    gcp => globalcell(:, 2, :)
                    frx => forcedrecv(:, 2, :)
                case (jmax)
                    ibp => iblank(:, jl, :)
                    gcp => globalcell(:, jl, :)
                    frx => forcedrecv(:, jl, :)
                case (kmin)
                    ibp => iblank(:, :, 2)
                    gcp => globalcell(:, :, 2)
                    frx => forcedrecv(:, :, 2)
                case (kmax)
                    ibp => iblank(:, :, kl)
                    gcp => globalcell(:, :, kl)
                    frx => forcedrecv(:, :, kl)
                end select
 
                ! -------------------------------------------------
                ! Step 1: Set the (haloed) cell iBlanks directly from
                ! the volume iBlanks
                ! -------------------------------------------------
                jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                ! Only do the slit elimination if we actually care about
                ! this surface for the zipper
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    ! -------------------------------------------------
                    ! Step 2: Slit elimination
                    ! -------------------------------------------------
 
                    ! Now we loop back through the cells again. For
                    ! interpolated cells with iblank=-1 we see if it satifies
                    ! the criteria above. If so we flag it wil "toFlip" =
                    ! 1. Note that we can't set a particular iblank directly
                    ! since that could cause an "avalance" effect with the
                    ! later iterations using the updated iblank from a previous
                    ! iteration.
                    allocate (toflip(ibeg:iend, jbeg:jend))
                    toflip = 0
                    do j = jbeg, jend
                        do i = ibeg, iend
 
                            ! We *might* add it if the interpolated cell is
                            ! touching two real cell on opposite sides.
 
                            if (bcdata(mm)%iBlank(i, j) == -1 .and. validcell(i, j)) then
 
                                ! Reset the side flag
                                side = .false.
 
                                if (validcell(i - 1, j) .and. bcdata(mm)%iBlank(i - 1, j) == 1) then
                                    side(1) = .true.
                                end if
 
                                if (validcell(i + 1, j) .and. bcdata(mm)%iBlank(i + 1, j) == 1) then
                                    side(2) = .true.
                                end if
 
                                if (validcell(i, j - 1) .and. bcdata(mm)%iBlank(i, j - 1) == 1) then
                                    side(3) = .true.
                                end if
 
                                if (validcell(i, j + 1) .and. bcdata(mm)%iBlank(i, j + 1) == 1) then
                                    side(4) = .true.
                                end if
 
                                if ((side(1) .and. side(2)) .or. (side(3) .and. side(4))) then
                                    toflip(i, j) = 1
                                end if
                            end if
                        end do
                    end do
 
                    ! Now just set the cell surface iblank to 1 if we
                    ! determined above we need to flip  the cell
 
                    do j = jbeg, jend
                        do i = ibeg, iend
                            if (toflip(i, j) == 1) then
                                bcdata(mm)%iBlank(i, j) = 1
                            end if
                        end do
                    end do
 
                    deallocate (toflip)
                end if faminclude
            end do bocoloop
        end do domainloop
 
    contains
        function validcell(i, j)
            implicit none
            integer(kind=intType), intent(in) :: i, j
            logical :: validcell
 
            ! for our purposes here, a valid cell is one that:
            ! 1. Is not a boundary halo. ie has globalCell >= 0
            ! 2. It is not a force receiver.
 
            validcell = .false.
            if (gcp(i + 1, j + 1) >= 0 .and. frx(i + 1, j + 1) == 0) then
                validcell = .true.
            end if
        end function validcell
    end subroutine slitelimination
 
    !
    !      surfaceDeviation computes an approximation of the maximum
    !      deviation a surface could be as compared to an underlying "exact"
    !      surface. The purpose is to compute an adaptive "near wall distance"
    !      value that can be used to determine if a point is "close" to a
    !     * wall.
    !
 
    subroutine surfacedeviation(famList, level, sps)
 
        use constants
        use blockpointers, only: bcdata, x, nbocos, ndom, bctype, il, jl, kl, bcfaceid
        use utils, only: setpointers, mynorm2, setbcpointers
        use sorting, only: faminlist
        use bcpointers, only: xx
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in), dimension(:) :: famList
        integer(kind=intType), intent(in) :: level, sps
 
        ! Local Variables
        integer(kind=intType) :: i, j, k, ii, jj, kk, nn, iBeg, iEnd, jBeg, jEnd, mm
        real(kind=realtype) :: deviation
 
        ! Loop over blocks
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            bocoloop: do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    call setbcpointers(mm, .true.)
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    ! The procedure goes in 2 passes. The first pass checks all
                    ! the i-direction edges, and the second all the j-direction
                    ! edges. For every edge, we estimate the max deviation
                    ! along that edge and then the surface will use the maximum
                    ! deviation from each of the 4 edges. Ie we scatter the
                    ! edge deviation to the two cells next to it. We only do
                    ! the real cells here. Boundary halos get -one set (below)
                    ! and then actual compute halos are set with an exchange.
 
                    bcdata(mm)%delta = -one
 
                    ! ------------------
                    ! Check the i-edges
                    ! ------------------
                    do j = jbeg, jend       ! <------- Node loop
                        do i = ibeg + 1, iend  ! <------- Face Loop
 
                            ! We will creating a local cubic approximation of the
                            ! local edge. This will use node i-2, i-1, i, and
                            ! i+1. However, due to the pointer offset, these are
                            ! all shifted by 1 to get: i-1, i, i+1, i+2
 
                            deviation = checkdeviation(xx(i - 1, j, :), xx(i, j, :), xx(i + 1, j, :), &
                                                       xx(i + 2, j, :))
 
                            ! Cell to the bottom:
                            if (j - 1 >= jbeg + 1) then
                                bcdata(mm)%delta(i, j - 1) = max(bcdata(mm)%delta(i, j - 1), deviation)
                            end if
 
                            ! Cell to the top:
                            if (j + 1 <= jend) then
                                bcdata(mm)%delta(i, j + 1) = max(bcdata(mm)%delta(i, j + 1), deviation)
                            end if
                        end do
                    end do
 
                    ! -----------------
                    ! Check the j-edges
                    ! -----------------
                    do j = jbeg + 1, jend   ! <------- Face loop
                        do i = ibeg, iend  ! <------- Node Loop
 
                            ! We will creating a local cubic approximation of the
                            ! local edge. This will use node j-2, j-1, j, and
                            ! j+1. However, due to the pointer offset, these are
                            ! all shifted by 1 to get: j-1, j, j+1, j+2
 
                            deviation = checkdeviation(xx(i, j - 1, :), xx(i, j, :), xx(i, j + 1, :), &
                                                       xx(i, j + 2, :))
 
                            ! Cell to the left:
                            if (i - 1 >= ibeg + 1) then
                                bcdata(mm)%delta(i - 1, j) = max(bcdata(mm)%delta(i - 1, j), deviation)
                            end if
 
                            ! Cell to the right:
                            if (i + 1 <= iend) then
                                bcdata(mm)%delta(i + 1, j) = max(bcdata(mm)%delta(i + 1, j), deviation)
                            end if
 
                        end do
                    end do
                end if faminclude
            end do bocoloop
        end do
 
        ! Exchange so that halos get correct values set as well. THIS IS BROKEN FIX IT!
        !call exchangeSurfaceDelta(level, sps, commPatternCell_1st, internalCell_1st)
 
        ! Now make one pass back and compute a delta for the nodes. Of
        ! course, this technically makes no sense: The nodes should
        ! exactly
        ! using this as a surrogate for what is near a surface, it make a
        ! little sense. Essentially we go through the nodes, and take the
        ! max deviation from the cells surrpounding it.
 
        ! Loop over blocks
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            bocoloop2: do mm = 1, nbocos
                faminclude2: if (faminlist(bcdata(mm)%famID, famlist)) 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
 
                            ! Since we are taking the max and the boundary halos
                            ! have a value of -one it's ok to blindy just take
                            ! the max from each of the 4 cells surrounding each node.
 
                            bcdata(mm)%deltaNode(i, j) = max( &
                                                         bcdata(mm)%delta(i, j), &
                                                         bcdata(mm)%delta(i + 1, j), &
                                                         bcdata(mm)%delta(i, j + 1), &
                                                         bcdata(mm)%delta(i + 1, j + 1))
                        end do
                    end do
                end if faminclude2
            end do bocoloop2
        end do
 
    end subroutine surfacedeviation
 
    function checkdeviation(P0, P1, P2, P3)
 
        ! Find the maximum deviation between a local cubic approximation
        ! formed by nodes P0, P1, P2 and P3, with the linear approximation
        ! formed by nodes P1 and P2.
 
        ! See this article for the implementation.
        ! https://en.wikipedia.org/wiki/Centripetal_Catmull-Rom_spline
 
        use constants
        use utils, only: mynorm2
        implicit none
 
        ! Input Parameters
        real(kind=realtype), intent(in), dimension(3) :: p0, p1, p2, p3
 
        ! Function value
        real(kind=realtype) :: checkdeviation
 
        ! Working Parameters
        real(kind=realtype) :: t0, t1, t2, t3
        real(kind=realtype), dimension(3) :: a1, a2, a3, b1, b2
        real(kind=realtype), parameter :: alpha = half
        integer(kind=intType), parameter :: n = 20
        integer(kind=intType) :: i
        real(kind=realtype) :: t, p(3), q(3), s
 
        t0 = zero
        t1 = t0 + mynorm2(p1 - p0)**alpha
        t2 = t1 + mynorm2(p2 - p1)**alpha
        t3 = t2 + mynorm2(p3 - p2)**alpha
 
        ! Normalize
        t1 = t1 / t3
        t2 = t2 / t3
        t3 = one
 
        ! Loop over the number of points to check. We need to go between t2
        ! and t3. No need to check the first and last since the devaition
        ! there is zero by construction.
        checkdeviation = zero
 
        do i = 1, n
            s = (i - one) / (n - one)
            t = (one - s) * t1 + s * t2
 
            ! Spline pt
            a3 = (t3 - t) / (t3 - t2) * p2 + (t - t2) / (t3 - t2) * p3
            a2 = (t2 - t) / (t2 - t1) * p1 + (t - t1) / (t2 - t1) * p2
            a1 = (t1 - t) / (t1 - t0) * p0 + (t - t0) / (t1 - t0) * p1
 
            b2 = (t3 - t) / (t3 - t1) * a2 + (t - t1) / (t3 - t1) * a3
            b1 = (t2 - t) / (t2 - t0) * a1 + (t - t0) / (t2 - t0) * a2
 
            p = (t2 - t) / (t2 - t1) * b1 + (t - t1) / (t2 - t1) * b2
 
            ! Now project the cubic point onto the line to get point Q
 
            q = p1 + dot_product(p - p1, p2 - p1) / dot_product(p2 - p1, p2 - p1) * (p2 - p1)
 
            ! Just get the distance between the two points.
            checkdeviation = max(checkdeviation, mynorm2(q - p))
 
        end do
 
    end function checkdeviation
 
    subroutine writewalls(famList)
 
        !use oversetData
        use constants
        use blockpointers
        use utils, only: setpointers, setbcpointers
        use bcpointers, only: xx
        use sorting, only: faminlist
        use communication, only: myid, adflow_comm_world, nproc
        use utils, only: echk
        use commonformats, only: sci12
        implicit none
        integer(kind=intType), intent(in), dimension(:) :: famList
        character(80) :: fileName, zoneName
        integer(kind=intType) :: i, j, nn, iDom, iBeg, iEnd, jBeg, jEnd, mm, iDim
        integer(kind=intType) :: nNode, nCell, tag, ierr, ii, iProc, nLocalBoco, tmp(5)
        real(kind=realtype), dimension(:), allocatable :: xbuffer
        integer(kind=intType), dimension(:), allocatable :: iblankBuffer, bocosPerProc
        integer(kind=intType), dimension(:, :, :), allocatable :: faceInfo
        integer, dimension(mpi_status_size) :: mpiStatus
 
        ! Write a gathered surface tecplot file.
        if (myid == 0) then
            write (filename, "(a)") "zipper_wall.dat"
 
            open (unit=101, file=trim(filename), form='formatted')
            write (101, *) 'TITLE = "zipper walls"'
            write (101, *) 'Variables = "X", "Y", "Z", "CellIBlank"'
        end if
 
        ! Before we start, the root processor needs to know how many
        ! receives we can expect from each processor.
 
        nlocalboco = 0
        do nn = 1, ndom
            call setpointers(nn, 1, 1)
            do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
                    nlocalboco = nlocalboco + 1
                end if faminclude
            end do
        end do
 
        if (myid == 0) then
            allocate (bocosperproc(0:nproc - 1))
        end if
 
        call mpi_gather(nlocalboco, 1, adflow_integer, bocosperproc, 1, &
                        adflow_integer, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
        ! Now setup the info array on the root proc:
        if (myid == 0) then
            allocate (faceinfo(5, maxval(bocosperproc), 0:nproc - 1))
        end if
 
        if (myid /= 0) then
            tag = 0
            do nn = 1, ndom
                call setpointers(nn, 1, 1)
                do mm = 1, nbocos
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                    faminclude2: if (faminlist(bcdata(mm)%famID, famlist)) then
                        tag = tag + 1
                        tmp = (/ibeg, iend, jbeg, jend, nbkglobal/)
                        call mpi_send(tmp, 5, &
                                      adflow_integer, 0, tag, adflow_comm_world, ierr)
                        call echk(ierr, __file__, __line__)
                    end if faminclude2
                end do
            end do
        else
            ! Receive the size info:
            do iproc = 1, nproc - 1
                do tag = 1, bocosperproc(iproc)
                    call mpi_recv(tmp, 5, adflow_integer, iproc, tag, &
                                  adflow_comm_world, mpistatus, ierr)
                    call echk(ierr, __file__, __line__)
                    faceinfo(:, tag, iproc) = tmp
                end do
            end do
        end if
 
        ! Need a barrier between the two sets of the comms just in case.
        call mpi_barrier(adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        tag = 0
        do nn = 1, ndom
            call setpointers(nn, 1, 1)
            do mm = 1, nbocos
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                faminclude3: if (faminlist(bcdata(mm)%famID, famlist)) then
                    call setbcpointers(mm, .true.)
 
                    nnode = (iend - ibeg + 1) * (jend - jbeg + 1)
                    ncell = (iend - ibeg + 1) * (jend - jbeg + 1)
                    allocate (iblankbuffer(ncell), xbuffer(3 * nnode))
                    ii = 0
                    do j = jbeg + 1, jend
                        do i = ibeg + 1, iend
                            ii = ii + 1
                            iblankbuffer(ii) = bcdata(mm)%iBlank(i, j)
                        end do
                    end do
 
                    ii = 0
                    do idim = 1, 3
                        do j = jbeg, jend
                            do i = ibeg, iend
                                ii = ii + 1
                                xbuffer(ii) = xx(i + 1, j + 1, idim)
                            end do
                        end do
                    end do
 
                    if (myid == 0) then
                        ! We can write it directly:
                        call writezone(ibeg, iend, jbeg, jend, nbkglobal, xbuffer, iblankbuffer)
                    else
 
                        tag = tag + 1
                        call mpi_send(iblankbuffer, ncell, adflow_integer, 0, tag, &
                                      adflow_comm_world, ierr)
                        call echk(ierr, __file__, __line__)
 
                        tag = tag + 1
                        call mpi_send(xbuffer, nnode * 3, adflow_real, 0, tag, &
                                      adflow_comm_world, ierr)
                        call echk(ierr, __file__, __line__)
 
                    end if
                    deallocate (iblankbuffer, xbuffer)
                end if faminclude3
            end do
        end do
 
        ! Complete the receives and writes on the root proc:
        if (myid == 0) then
            ! Receive the nodes and iblank info
            do iproc = 1, nproc - 1
                do tag = 1, bocosperproc(iproc)
 
                    ibeg = faceinfo(1, tag, iproc)
                    iend = faceinfo(2, tag, iproc)
 
                    jbeg = faceinfo(3, tag, iproc)
                    jend = faceinfo(4, tag, iproc)
 
                    nbkglobal = faceinfo(5, tag, iproc)
 
                    nnode = (iend - ibeg + 1) * (jend - jbeg + 1)
                    ncell = (iend - ibeg + 1) * (jend - jbeg + 1)
 
                    allocate (iblankbuffer(ncell), xbuffer(3 * nnode))
 
                    call mpi_recv(iblankbuffer, ncell, adflow_integer, iproc, (2 * tag - 1), &
                                  adflow_comm_world, mpistatus, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call mpi_recv(xbuffer, nnode * 3, adflow_real, iproc, 2 * tag, &
                                  adflow_comm_world, mpistatus, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call writezone(ibeg, iend, jbeg, jend, nbkglobal, xbuffer, iblankbuffer)
 
                    deallocate (iblankbuffer, xbuffer)
 
                end do
            end do
            deallocate (faceinfo, bocosperproc)
            close (101)
        end if
 
    contains
 
        subroutine writezone(iBeg, iEnd, jBeg, jEnd, nBkGlobal, xx, iblank)
 
            implicit none
            ! Input
            integer(kind=intType), intent(in) :: iBeg, iEnd, jBeg, jEnd, nBkGlobal
            real(kind=realtype), intent(in), dimension(:) :: xx
            integer(kind=intType), intent(in), dimension(:) :: iblank
 
            character(80) :: zoneName
            integer(kind=intType) :: iDim
            character(len=maxStringLen) :: zoneFormat
 
            write (zonename, "(a,I5.5)") "Zone_", nbkglobal
            zoneformat = "(3(A), I5, A, I5)"
            write (101, zoneformat) 'ZONE T=', trim(zonename), " I=", iend - ibeg + 1, " J=", jend - jbeg + 1
            write (101, *) "DATAPACKING=BLOCK, VARLOCATION=([1,2,3]=NODAL, [4]=CELLCENTERED)"
 
            ! The 3 is for the three coordinate directions
            nnode = (iend - ibeg + 1) * (jend - jbeg + 1)
            ncell = (iend - ibeg) * (jend - jbeg)
 
            do i = 1, 3 * nnode
                write (101, sci12) xx(i)
            end do
 
            do i = 1, ncell
                write (101, *) iblank(i)
            end do
        end subroutine writezone
    end subroutine writewalls
 
    subroutine bowtieandisolationelimination(famList, level, sps)
 
        use constants
        use blockpointers
        use communication
        use utils, only: setpointers
        use oversetcommutilities, only: exchangesurfaceiblanks
        use sorting, only: faminlist
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in), dimension(:) :: famList
        integer(kind=intType), intent(in) :: level, sps
 
        ! Local variables
        integer(kind=intType) :: mm, nn, i, j, k, e, iBeg, iEnd, jBeg, jEnd
        logical :: side(4)
 
        integer(kind=intType), dimension(:, :), pointer :: ibp, gcp
        integer(kind=intType), dimension(:, :), allocatable :: toFlip, nE, nC
 
        ! This routine initializes the surface cell iblank based on the
        ! volume iblank. It is not a straight copy since we a little
        ! preprocessing to eliminate a few particularly nasty cases.
        ! Three analysis are performed:
        ! 1. Bow-tie elimination
        ! 2. Single cell elmination
 
        bowtieloop: do e = 0, 2
            domainloop1: do nn = 1, ndom
                call setpointers(nn, level, sps)
 
                bocoloop1: do mm = 1, nbocos
                    faminclude1: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                        select case (bcfaceid(mm))
                        case (imin)
                            ibp => iblank(2, :, :)
                            gcp => globalcell(2, :, :)
                        case (imax)
                            ibp => iblank(il, :, :)
                            gcp => globalcell(il, :, :)
                        case (jmin)
                            ibp => iblank(:, 2, :)
                            gcp => globalcell(:, 2, :)
                        case (jmax)
                            ibp => iblank(:, jl, :)
                            gcp => globalcell(:, jl, :)
                        case (kmin)
                            ibp => iblank(:, :, 2)
                            gcp => globalcell(:, :, 2)
                        case (kmax)
                            ibp => iblank(:, :, kl)
                            gcp => globalcell(:, :, kl)
                        end select
 
                        ! -------------------------------------------------
                        ! Step 2: Bow-tie elimination: Elimiate cells
                        ! that touch only at a corner.
                        ! -------------------------------------------------
 
                        ! Make bounds a little easier to read. Owned cells only
                        ! from now on.
                        jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                        ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                        ! Allocate two tmporary auxilary arrays 'eN'->
                        ! edgeNeighbours and 'cN'-> cornerNeighbous. For every
                        ! comupte cell determine the number of compute
                        ! neighbours connected along edges and at corners
                        allocate (ne(ibeg:iend, jbeg:jend), nc(ibeg:iend, jbeg:jend))!, &
 
                        call findbowties()
 
                        do j = jbeg, jend
                            do i = ibeg, iend
                                if (bcdata(mm)%iBlank(i, j) > 0 .and. nc(i, j) >= 1 .and. ne(i, j) <= e) then
                                    bcdata(mm)%iBlank(i, j) = 0
                                end if
                            end do
                        end do
 
                        deallocate (nc, ne)
                    end if faminclude1
                end do bocoloop1
            end do domainloop1
 
            ! Since we potentially changed iBlanks, we need to updated by
            ! performing an exchange.
            call exchangesurfaceiblanks(famlist, level, sps, commpatterncell_2nd, internalcell_2nd)
 
        end do bowtieloop
 
        domainloop2: do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            bocoloop2: do mm = 1, nbocos
                faminclude2: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        gcp => globalcell(2, :, :)
                    case (imax)
                        gcp => globalcell(il, :, :)
                    case (jmin)
                        gcp => globalcell(:, 2, :)
                    case (jmax)
                        gcp => globalcell(:, jl, :)
                    case (kmin)
                        gcp => globalcell(:, :, 2)
                    case (kmax)
                        gcp => globalcell(:, :, kl)
                    end select
 
                    ! Make bounds a little easier to read. Owned cells only
                    ! from now on.
                    jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                    ! -------------------------------------------------
                    ! Step 3: Single-cell elimination: Elimiate cells
                    ! that do not touch any other cells.
                    ! -------------------------------------------------
 
                    allocate (ne(ibeg:iend, jbeg:jend), nc(ibeg:iend, jbeg:jend))
 
                    call setneighbourcounts()
 
                    ! This is easy, if a compute cell is stil around with no
                    ! neighbours, kill it
                    do j = jbeg, jend
                        do i = ibeg, iend
                            if (bcdata(mm)%iBlank(i, j) == 1 .and. ne(i, j) == 0 .and. nc(i, j) == 0) then
                                bcdata(mm)%iBlank(i, j) = 0
                            end if
                        end do
                    end do
 
                    deallocate (ne, nc)
 
                end if faminclude2
            end do bocoloop2
        end do domainloop2
 
        ! Again, since we potentially changed iBlanks, we need to updated by
        ! performing an exchange.
        call exchangesurfaceiblanks(famlist, level, sps, commpatterncell_2nd, internalcell_2nd)
    contains
        subroutine findbowties
 
            implicit none
            ! For every compute determine the number of compute neighbours
            ! connected along edges (nE) and the number of bow-ties (in nC)
 
            integer(kind=intType) :: i, j, e(4)
 
            ne = 0
            nc = 0
 
            do j = jbeg, jend
                do i = ibeg, iend
                    if (bcdata(mm)%iBlank(i, j) >= 0) then
                        e = 0
 
                        !    |  e3  |
                        !  --c4-----c3-
                        !    |      |
                        ! e4 |   x  | e2
                        !    |      |
                        ! -- c1-----c2-
                        !    |  e1  |
 
                        ! Set the status of each of the 4 edges:
 
                        if (bcdata(mm)%iBlank(i, j - 1) == 1) &
                            e(1) = 1
 
                        if (bcdata(mm)%iBlank(i + 1, j) == 1) &
                            e(2) = 1
 
                        if (bcdata(mm)%iBlank(i, j + 1) == 1) &
                            e(3) = 1
 
                        if (bcdata(mm)%iBlank(i - 1, j) == 1) &
                            e(4) = 1
 
                        ! Check the 4 corner neighbours for bow-tie status
                        if (bcdata(mm)%iBlank(i - 1, j - 1) == 1 .and. e(4) == 0 .and. e(1) == 0) &
                            nc(i, j) = nc(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i + 1, j - 1) == 1 .and. e(1) == 0 .and. e(2) == 0) &
                            nc(i, j) = nc(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i + 1, j + 1) == 1 .and. e(2) == 0 .and. e(3) == 0) &
                            nc(i, j) = nc(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i - 1, j + 1) == 1 .and. e(3) == 0 .and. e(4) == 0) &
                            nc(i, j) = nc(i, j) + 1
 
                        ne(i, j) = sum(e)
                    end if
                end do
            end do
        end subroutine findbowties
 
        subroutine setneighbourcounts
 
            implicit none
            ! For every comute determine the number of compute neighbours
            ! connected along edges and at corners
 
            integer(kind=intType) :: i, j
 
            ne = 0
            nc = 0
 
            do j = jbeg, jend
                do i = ibeg, iend
                    if (bcdata(mm)%iBlank(i, j) == 1) then
 
                        !    |  e3  |
                        !  --c4-----c3-
                        !    |      |
                        ! e4 |   x  | e2
                        !    |      |
                        ! -- c1-----c2-
                        !    |  e1  |
 
                        ! Set the status of each of the 4 edges:
 
                        if (bcdata(mm)%iBlank(i, j - 1) == 1) &
                            ne(i, j) = ne(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i + 1, j) == 1) &
                            ne(i, j) = ne(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i, j + 1) == 1) &
                            ne(i, j) = ne(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i - 1, j) == 1) &
                            ne(i, j) = ne(i, j) + 1
 
                        ! Check the 4 corner neighbours for compute neighbour
                        if (bcdata(mm)%iBlank(i - 1, j - 1) == 1) &
                            nc(i, j) = nc(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i + 1, j - 1) == 1) &
                            nc(i, j) = nc(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i + 1, j + 1) == 1) &
                            nc(i, j) = nc(i, j) + 1
 
                        if (bcdata(mm)%iBlank(i - 1, j + 1) == 1) &
                            nc(i, j) = nc(i, j) + 1
                    end if
                end do
            end do
        end subroutine setneighbourcounts
    end subroutine bowtieandisolationelimination
end module zippermesh