makeBoundaryStrings.F90

Go to the documentation of this file.

module gapboundaries
contains
    subroutine makegapboundarystrings(zipperFamList, level, sps, master)
 
        use constants
        use adtbuild, only: buildserialquad
        use blockpointers, only: x, globalcell, globalnode, bcdata, nbocos, &
                                 il, jl, kl, ndom, righthanded, bcfaceid, bctype
        use communication, only: adflow_comm_world, myid, nproc
        use oversetdata, only: oversetstring, oversetwall, nclusters, clusters, cumdomproc
        use stringops
        use utils, only: setpointers, echk, mynorm2, cross_prod
        use oversetpackingroutines, only: getwallsize
        use sorting, only: faminlist
        implicit none
 
        ! Input Params
        integer(kind=intType), intent(in), dimension(:) :: zipperFamLIst
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working
        integer(kind=intType) :: i, j, k, nn, mm, ii, jj, kk, c, e, idx
        integer(kind=intType) :: i1, i2, j1, j2, iBeg, iEnd, jBeg, jEnd
        integer(kind=intType) :: i3, i4, j3, j4
        integer(kind=intType) :: iStart, iSize, ierr, iProc, firstElem, curElem
        integer(kind=intType) :: below, above, left, right, nNodes, nElems
        integer(kind=intType) :: patchNodeCounter, nZipped, gc
        integer(kind=intType), dimension(:), allocatable :: nElemsProc, nNodesProc
        integer(kind=intType), dimension(:, :), pointer :: gcp
        real(kind=realtype), dimension(:, :, :), pointer :: xx
        real(kind=realtype), dimension(3) :: s1, s2, s3, s4, v1, v2, v3, v4, x0
        real(kind=realtype) :: fact, timea, minnorm
 
        real(kind=realtype), dimension(:, :, :), allocatable :: patchnormals
        real(kind=realtype), dimension(:, :), allocatable :: patchh
        integer(kind=intType), dimension(:), allocatable :: epc, surfaceSeeds, inverse
        logical, dimension(:), allocatable :: badString
        type(oversetstring), dimension(:), allocatable, target :: localStrings
        type(oversetstring), dimension(:), allocatable, target :: globalStrings
        type(oversetstring) :: master, pocketMaster
        type(oversetstring), pointer :: stringsLL, str
        type(oversetstring), dimension(:), allocatable, target :: strings
        integer(kind=intType) :: nFullStrings, nUnique, famID
        logical :: regularOrdering
        integer mpiStatus(MPI_STATUS_SIZE)
 
        ! Wall search related
        integer(kind=intType) :: ncells
        type(oversetwall), dimension(:), allocatable, target :: walls
        type(oversetwall), target :: fullWall
        character(80) :: fileName
 
        ! Set small number to avoid division by zero when computing normal vectors
        minnorm = 1.0e-14
 
        ! Loop over the wall faces counting up the edges that stradle a
        ! compute cell and a blanked (or interpolated) cell.
 
        allocate (epc(nclusters)) ! epc = elementsPerCluster
        epc = 0
 
        ! Get a (very) large overestimate of the total number of edges in a
        ! cluster: as twice the  number of nodes.
        domainloop: do nn = 1, ndom
            call setpointers(nn, level, sps)
            call getwallsize(zipperfamlist, nnodes, nelems, .false.)
            c = clusters(cumdomproc(myid) + nn)
            epc(c) = epc(c) + 2 * nnodes
        end do domainloop
 
        ! Allocate the space we need in the local strings.
        allocate (localstrings(nclusters))
        do c = 1, nclusters
            call nullifystring(localstrings(c))
        end do
 
        do c = 1, nclusters
            allocate ( &
                localstrings(c)%conn(2, epc(c)), localstrings(c)%nodeData(10, 2 * epc(c)), &
                localstrings(c)%intNodeData(3, 2 * epc(c)))
            localstrings(c)%nodeData = zero
            localstrings(c)%nNodes = 0
            localstrings(c)%nElems = 0
 
            ! Assign string pointers immediately after allocation
            call setstringpointers(localstrings(c))
 
        end do
        deallocate (epc)
        ! And now loop back through the walls and add in the
        ! elems/nodes/normals/indices for each edge.
 
        ! Reset the elems per cluster. We will count up the actual number
        ! now.
 
        patchnodecounter = 0
        domainloop2: do nn = 1, ndom
            call setpointers(nn, level, sps)
            ! The current cluster is 'c'
            c = clusters(cumdomproc(myid) + nn)
 
            bocoloop: do mm = 1, nbocos
                famid = bcdata(mm)%famID
                if (faminlist(famid, zipperfamlist)) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        xx => x(1, :, :, :)
                        gcp => globalcell(2, :, :)
                        fact = one
                        regularordering = .true.
                    case (imax)
                        xx => x(il, :, :, :)
                        gcp => globalcell(il, :, :)
                        fact = -one
                        regularordering = .false.
                    case (jmin)
                        xx => x(:, 1, :, :)
                        gcp => globalcell(:, 2, :)
                        fact = -one
                        regularordering = .false.
                    case (jmax)
                        xx => x(:, jl, :, :)
                        gcp => globalcell(:, jl, :)
                        fact = one
                        regularordering = .true.
                    case (kmin)
                        xx => x(:, :, 1, :)
                        gcp => globalcell(:, :, 2)
                        fact = one
                        regularordering = .true.
                    case (kmax)
                        xx => x(:, :, kl, :)
                        gcp => globalcell(:, :, kl)
                        fact = -one
                        regularordering = .false.
                    end select
 
                    ! Need to reverse once more for a left-handed block
                    if (.not. righthanded) then
                        fact = -fact
                        regularordering = .not. (regularordering)
                    end if
 
                    ! Before we go through and find the actual elems,
                    ! precompute the patch numbering and node-based averaged
                    ! unit normals
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    allocate (patchnormals(3, ibeg:iend, jbeg:jend), &
                              patchh(ibeg:iend, jbeg:jend))
 
                    do j = jbeg, jend
                        do i = ibeg, iend
                            patchnodecounter = patchnodecounter + 1
                            x0 = xx(i + 1, j + 1, :)
 
                            ! Normalized normals for each surrounding face.
                            v1 = xx(i + 2, j + 1, :) - x0
                            v2 = xx(i + 1, j + 2, :) - x0
                            v3 = xx(i, j + 1, :) - x0
                            v4 = xx(i + 1, j, :) - x0
 
                            call cross_prod(v1, v2, s1)
                            call cross_prod(v2, v3, s2)
                            call cross_prod(v3, v4, s3)
                            call cross_prod(v4, v1, s4)
 
                            ! When we have an 0-grid node, two of the v vectors will be the same.
                            ! Therefore, one of the s vectors will be zero. So we define minNorm
                            ! to avoid a division by zero. This will not affect the averaged normal
                            ! vector, since we will normalize it anyway in the end.
 
                            s1 = s1 / max(minnorm, mynorm2(s1))
                            s2 = s2 / max(minnorm, mynorm2(s2))
                            s3 = s3 / max(minnorm, mynorm2(s3))
                            s4 = s4 / max(minnorm, mynorm2(s4))
 
                            ! Average and do final normalization including
                            ! correcting for inward normals.
                            s1 = fourth * (s1 + s2 + s3 + s4)
                            patchnormals(:, i, j) = s1 / mynorm2(s1) * fact
 
                            ! Get the maximum edge length for this node. Use the
                            ! 4 diagonal nodes:
                            v1 = xx(i + 2, j + 2, :) - x0
                            v2 = xx(i, j + 2, :) - x0
                            v3 = xx(i, j, :) - x0
                            v4 = xx(i + 2, j, :) - x0
 
                            patchh(i, j) = max(mynorm2(v1), mynorm2(v2), mynorm2(v3), mynorm2(v4))
 
                        end do
                    end do
 
                    ! ------------------
                    ! Check the i-edges
                    ! ------------------
                    do j = jbeg, jend       ! <------- Node loop
                        do i = ibeg + 1, iend  ! <------- Face Loop
                            if (gcp(i + 1, j + 1) >= 0 .and. gcp(i + 1, j + 2) >= 0) then
                                below = max(bcdata(mm)%iBlank(i, j), 0)
                                above = max(bcdata(mm)%iBlank(i, j + 1), 0)
 
                                if ((below == 0 .and. above == 1) .or. (below == 1 .and. above == 0)) then
                                    localstrings(c)%nNodes = localstrings(c)%nNodes + 2
                                    localstrings(c)%nElems = localstrings(c)%nElems + 1
                                    e = localstrings(c)%nElems
 
                                    ! Make sure the real cell is on the LEFT
                                    ! of the directed edge
 
                                    if (below == 0) then
                                        i1 = i - 1; j1 = j
                                        i2 = i; j2 = j
 
                                        i3 = i1; j3 = j + 1
                                        i4 = i2; j4 = j + 1
                                    else
                                        i1 = i; j1 = j
                                        i2 = i - 1; j2 = j
 
                                        i3 = i1; j3 = j - 1
                                        i4 = i2; j4 = j - 1
                                    end if
 
                                    ! Don't forget pointer offset for xx
                                    if (regularordering) then
                                        localstrings(c)%nodeData(1:3, 2 * e - 1) = xx(i1 + 1, j1 + 1, :)
                                        localstrings(c)%nodeData(1:3, 2 * e) = xx(i2 + 1, j2 + 1, :)
 
                                        ! Global index of node on reduced global surface.
                                        localstrings(c)%intNodeData(1, 2 * e - 1) = bcdata(mm)%surfIndex(i1, j1)
                                        localstrings(c)%intNodeData(1, 2 * e) = bcdata(mm)%surfIndex(i2, j2)
                                    else
                                        localstrings(c)%nodeData(1:3, 2 * e) = xx(i1 + 1, j1 + 1, :)
                                        localstrings(c)%nodeData(1:3, 2 * e - 1) = xx(i2 + 1, j2 + 1, :)
 
                                        localstrings(c)%intNodeData(1, 2 * e) = bcdata(mm)%surfIndex(i1, j1)
                                        localstrings(c)%intNodeData(1, 2 * e - 1) = bcdata(mm)%surfIndex(i2, j2)
 
                                    end if
                                    v1 = xx(i1 + 1, j1 + 1, :) - xx(i3 + 1, j3 + 1, :)
                                    v1 = v1 / mynorm2(v1)
 
                                    v2 = xx(i2 + 1, j2 + 1, :) - xx(i4 + 1, j4 + 1, :)
                                    v2 = v2 / mynorm2(v2)
 
                                    ! Perpendicular vector
                                    localstrings(c)%nodeData(7:9, 2 * e - 1) = v1
                                    localstrings(c)%nodeData(7:9, 2 * e) = v2
 
                                    ! Averaged node normal
                                    localstrings(c)%nodeData(4:6, 2 * e - 1) = patchnormals(:, i1, j1)
                                    localstrings(c)%nodeData(4:6, 2 * e) = patchnormals(:, i2, j2)
 
                                    ! Surface deviation estimation
                                    localstrings(c)%nodeData(10, 2 * e - 1) = patchh(i1, j1)
                                    localstrings(c)%nodeData(10, 2 * e) = patchh(i2, j2)
 
                                    ! Cluster of the node
                                    localstrings(c)%intNodeData(2, 2 * e - 1) = c
                                    localstrings(c)%intNodeData(2, 2 * e) = c
 
                                    ! Family ID of node. Not implemented yet.
                                    localstrings(c)%intNodeData(3, 2 * e - 1) = famid
                                    localstrings(c)%intNodeData(3, 2 * e) = famid
 
                                    ! Connectivity
                                    localstrings(c)%conn(:, e) = (/2 * e - 1, 2 * e/)
                                end if
                            end if
                        end do
                    end do
 
                    ! -----------------
                    ! Check the j-edges
                    ! -----------------
                    do j = jbeg + 1, jend   ! <------- Face loop
                        do i = ibeg, iend  ! <------- Node Loop
                            if (gcp(i + 1, j + 1) >= 0 .and. gcp(i + 2, j + 1) >= 0) then
                                left = max(bcdata(mm)%iBlank(i, j), 0)
                                right = max(bcdata(mm)%iBlank(i + 1, j), 0)
 
                                if ((left == 0 .and. right == 1) .or. (left == 1 .and. right == 0)) then
                                    localstrings(c)%nNodes = localstrings(c)%nNodes + 2
                                    localstrings(c)%nElems = localstrings(c)%nElems + 1
 
                                    e = localstrings(c)%nElems
 
                                    ! Again, make sure the real cell is on the LEFT
                                    ! of the directed edge
                                    if (left == 0) then
                                        i1 = i; j1 = j
                                        i2 = i; j2 = j - 1
 
                                        i3 = i1 + 1; j3 = j1
                                        i4 = i2 + 1; j4 = j2
 
                                    else
                                        i1 = i; j1 = j - 1
                                        i2 = i; j2 = j
 
                                        i3 = i1 - 1; j3 = j1
                                        i4 = i2 - 1; j4 = j2
                                    end if
 
                                    ! Don't forget pointer offset xx
                                    if (regularordering) then
                                        localstrings(c)%nodeData(1:3, 2 * e - 1) = xx(i1 + 1, j1 + 1, :)
                                        localstrings(c)%nodeData(1:3, 2 * e) = xx(i2 + 1, j2 + 1, :)
 
                                        ! Index of global node
                                        localstrings(c)%intNodeData(1, 2 * e - 1) = bcdata(mm)%surfIndex(i1, j1)
                                        localstrings(c)%intNodeData(1, 2 * e) = bcdata(mm)%surfIndex(i2, j2)
 
                                    else
                                        localstrings(c)%nodeData(1:3, 2 * e) = xx(i1 + 1, j1 + 1, :)
                                        localstrings(c)%nodeData(1:3, 2 * e - 1) = xx(i2 + 1, j2 + 1, :)
 
                                        ! Index of global node
                                        localstrings(c)%intNodeData(1, 2 * e) = bcdata(mm)%surfIndex(i1, j1)
                                        localstrings(c)%intNodeData(1, 2 * e - 1) = bcdata(mm)%surfIndex(i2, j2)
 
                                    end if
 
                                    v1 = xx(i1 + 1, j1 + 1, :) - xx(i3 + 1, j3 + 1, :)
                                    v1 = v1 / mynorm2(v1)
 
                                    v2 = xx(i2 + 1, j2 + 1, :) - xx(i4 + 1, j4 + 1, :)
                                    v2 = v2 / mynorm2(v2)
 
                                    ! Perpendicular vector
                                    localstrings(c)%nodeData(7:9, 2 * e - 1) = v1
                                    localstrings(c)%nodeData(7:9, 2 * e) = v2
 
                                    ! Averaged node normal
                                    localstrings(c)%nodeData(4:6, 2 * e - 1) = patchnormals(:, i1, j1)
                                    localstrings(c)%nodeData(4:6, 2 * e) = patchnormals(:, i2, j2)
 
                                    ! Surface deviation estimation
                                    localstrings(c)%nodeData(10, 2 * e - 1) = patchh(i1, j1)
                                    localstrings(c)%nodeData(10, 2 * e) = patchh(i2, j2)
 
                                    ! Cluster of the node
                                    localstrings(c)%intNodeData(2, 2 * e - 1) = c
                                    localstrings(c)%intNodeData(2, 2 * e) = c
 
                                    ! Family ID of node. Not implemented yet.
                                    localstrings(c)%intNodeData(3, 2 * e - 1) = famid
                                    localstrings(c)%intNodeData(3, 2 * e) = famid
 
                                    ! Connectivity
                                    localstrings(c)%conn(:, e) = (/2 * e - 1, 2 * e/)
                                end if
                            end if
                        end do
                    end do
                    deallocate (patchnormals, patchh)
                end if
            end do bocoloop
        end do domainloop2
 
        ! Before we send the gap strings to the root proc, reduce them so
        ! the root proc has a little less work to do.
        do c = 1, nclusters
            call reducegapstring(localstrings(c))
        end do
 
        ! Allocate the global list of strings on the root proc
        if (myid == 0) then
            allocate (globalstrings(nclusters))
            do c = 1, nclusters
                call nullifystring(globalstrings(c))
            end do
        end if
 
        ! Next for each each cluster, gather to the root the gap boundary strings
 
        allocate (nnodesproc(0:nproc), nelemsproc(0:nproc))
 
        do c = 1, nclusters
            ! Now let the root processor know how many nodes/elements my
            ! processor will be sending:
            nelemsproc(0) = 0
            nnodesproc(0) = 0
 
            call mpi_gather(localstrings(c)%nElems, 1, adflow_integer, nelemsproc(1:nproc), 1, adflow_integer, 0, &
                            adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            call mpi_gather(localstrings(c)%nNodes, 1, adflow_integer, nnodesproc(1:nproc), 1, adflow_integer, 0, &
                            adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            if (myid == 0) then
 
                ! Before we can receive stuff, we need to determine the node
                ! off-sets such that the conn from the strings on each processor
                ! don't overlap.
 
                do i = 2, nproc
                    ! The 0 and 1st entry of the nEdgeProc and nNodeProc arrays are already correct:
                    nnodesproc(i) = nnodesproc(i) + nnodesproc(i - 1)
                    nelemsproc(i) = nelemsproc(i) + nelemsproc(i - 1)
                end do
 
                allocate (globalstrings(c)%nodeData(10, nnodesproc(nproc)), &
                          globalstrings(c)%intNodeData(3, nnodesproc(nproc)), &
                          globalstrings(c)%conn(2, nelemsproc(nproc)))
 
                ! Always set the pointers immediately after allocation
                call setstringpointers(globalstrings(c))
 
                ! Put proc 0's own nodes/normals/indices in the global list if we have any
                do i = 1, localstrings(c)%nNodes
                    globalstrings(c)%nodeData(:, i) = localstrings(c)%nodeData(:, i)
                    globalstrings(c)%intNodeData(:, i) = localstrings(c)%intNodeData(:, i)
                end do
 
                ! Put proc 0's own elements in the global list if we have any
                do i = 1, localstrings(c)%nElems
                    globalstrings(c)%conn(:, i) = localstrings(c)%conn(:, i)
                end do
 
                ! Set my total sizes
                globalstrings(c)%nNodes = nnodesproc(nproc)
                globalstrings(c)%nElems = nelemsproc(nproc)
 
                ! Now receive from each of the other procs.
                do iproc = 1, nproc - 1
                    ! Check if this proc actually has anything to send:
                    if ((nelemsproc(iproc + 1) - nelemsproc(iproc)) > 0) then
                        istart = nnodesproc(iproc) + 1
                        iend = nnodesproc(iproc + 1)
                        isize = iend - istart + 1
 
                        ! ----------- Node sized arrays -------------
                        call mpi_recv(globalstrings(c)%nodeData(:, istart:iend), isize * 10, &
                                      adflow_real, iproc, iproc, &
                                      adflow_comm_world, mpistatus, ierr)
                        call echk(ierr, __file__, __line__)
 
                        call mpi_recv(globalstrings(c)%intNodeData(:, istart:iend), isize * 3, &
                                      adflow_integer, iproc, iproc, &
                                      adflow_comm_world, mpistatus, ierr)
                        call echk(ierr, __file__, __line__)
 
                        ! ----------- Element sized arrays -------------
                        istart = nelemsproc(iproc) + 1
                        iend = nelemsproc(iproc + 1)
                        isize = iend - istart + 1
                        call mpi_recv(globalstrings(c)%conn(:, istart:iend), isize * 2, adflow_integer, iproc, iproc, &
                                      adflow_comm_world, mpistatus, ierr)
                        call echk(ierr, __file__, __line__)
 
                        ! Increment the conn we just received by the node offset:
                        do i = istart, iend
                            globalstrings(c)%conn(:, i) = globalstrings(c)%conn(:, i) + nnodesproc(iproc)
                        end do
                    end if
                end do
            else
                ! Not root proc so send my stuff if we have anything:
                if (localstrings(c)%nElems > 0) then
 
                    ! ----------- Node sized arrays -------------
                    call mpi_send(localstrings(c)%nodeData, 10 * localstrings(c)%nNodes, adflow_real, 0, myid, &
                                  adflow_comm_world, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call mpi_send(localstrings(c)%intNodeData, 3 * localstrings(c)%nNodes, adflow_integer, 0, myid, &
                                  adflow_comm_world, ierr)
                    call echk(ierr, __file__, __line__)
 
                    ! ----------- Element sized arrays -------------
                    call mpi_send(localstrings(c)%conn, 2 * localstrings(c)%nElems, adflow_integer, 0, myid, &
                                  adflow_comm_world, ierr)
                    call echk(ierr, __file__, __line__)
 
                end if
            end if
        end do
 
        ! Everyone is now done with the local strings
        do c = 1, nclusters
            call deallocatestring(localstrings(c))
        end do
        deallocate (localstrings)
 
        ! Before we perform serial code implementations, surface wall info
        ! need to be communicated to root proc (zero) too. This will be used
        ! later to identify zipper triangle containment search for identifying
        ! quad surface cell info for force integration. Search will be
        ! performed on dual surface cells.
 
        ! ---------- Begin wall data accumulation -------------
        allocate (walls(nclusters))
        ! Build primal quad walls ADT
        call buildclusterwalls(level, sps, .false., walls, zipperfamlist, size(zipperfamlist))
 
        ! Finally build up a "full wall" that is made up of all the cluster
        ! walls.
 
        nnodes = 0
        ncells = 0
        do i = 1, nclusters
            nnodes = nnodes + walls(i)%nNodes
            ncells = ncells + walls(i)%nCells
        end do
 
        allocate (fullwall%x(3, nnodes))
        allocate (fullwall%conn(4, ncells))
        allocate (fullwall%ind(nnodes))
        allocate (fullwall%indCell(ncells))
 
        nnodes = 0
        ncells = 0
        ii = 0
        do i = 1, nclusters
 
            ! Add in the nodes/elements from this cluster
 
            do j = 1, walls(i)%nNodes
                nnodes = nnodes + 1
                fullwall%x(:, nnodes) = walls(i)%x(:, j)
                fullwall%ind(nnodes) = walls(i)%ind(j)
            end do
 
            do j = 1, walls(i)%nCells
                ncells = ncells + 1
                fullwall%conn(:, ncells) = walls(i)%conn(:, j) + ii
                fullwall%indCell(ncells) = walls(i)%indCell(j)
            end do
 
            ! Increment the node offset
            ii = ii + walls(i)%nNodes
        end do
 
        ! Finish the setup of the full wall.
        fullwall%nCells = ncells
        fullwall%nNodes = nnodes
        call buildserialquad(ncells, nnodes, fullwall%x, fullwall%conn, fullwall%ADT)
 
        ! Now all the procs have fullWall info.  Note: this is overkill,
        ! since only proc 0 needs them for zipper triangle containment
        ! search.
 
        ! =================================================================
        !                   Serial code from here on out
        ! =================================================================
 
        if (myid == 0) then
            timea = mpi_wtime()
 
            ! First thing we do is reduce each of the global cluster gap
            ! strings
            do c = 1, nclusters
                call reducegapstring(globalstrings(c))
            end do
 
            ! Combine all global strings together into a masterString. First
            ! count up the sizes
            nelems = 0
            nnodes = 0
            do c = 1, nclusters
                nelems = nelems + globalstrings(c)%nElems
                nnodes = nnodes + globalstrings(c)%nNodes
            end do
 
            call nullifystring(master)
            master%nNodes = nnodes
            master%nElems = nelems
            allocate (master%nodeData(10, nnodes), master%conn(2, nelems), &
                      master%intNodeData(3, nnodes))
 
            ! Set the string pointers to the individual arrays
            call setstringpointers(master)
 
            nnodes = 0 ! This is our running counter for offseting nodes
            ii = 0
            jj = 0
 
            do c = 1, nclusters
                do i = 1, globalstrings(c)%nNodes
                    ii = ii + 1
                    master%nodeData(:, ii) = globalstrings(c)%nodeData(:, i)
                    master%intNodeData(:, ii) = globalstrings(c)%intNodeData(:, i)
                end do
 
                do i = 1, globalstrings(c)%nElems
                    jj = jj + 1
                    master%conn(:, jj) = globalstrings(c)%conn(:, i) + nnodes
                end do
                nnodes = ii
            end do
 
            ! Now the root is done with the global strings so deallocate that
            ! too.
            do c = 1, nclusters
                call deallocatestring(globalstrings(c))
            end do
            deallocate (globalstrings)
        end if
 
    end subroutine makegapboundarystrings
end module gapboundaries