stringOps.F90

Go to the documentation of this file.

module stringops
 
    ! Import oversetString becuase every routine uses this.
    use oversetdata, only: oversetstring
contains
 
    subroutine nullifystring(string)
 
        use constants
        implicit none
 
        type(oversetstring) :: string
 
        nullify (string%nodeData, &
                 string%x, &
                 string%norm, &
                 string%perpNorm, &
                 string%h, &
                 string%intNodeData, &
                 string%ind, &
                 string%cluster, &
                 string%family, &
                 string%conn, &
                 string%pNodes, &
                 string%pElems, &
                 string%cNodes, &
                 string%otherID, &
                 string%nte, &
                 string%subStr, &
                 string%elemUsed, &
                 string%XzipNodeUsed, &
                 string%tris, &
                 string%surfCellID)
 
    end subroutine nullifystring
 
    subroutine deallocatestring(string)
 
        use constants
        implicit none
 
        type(oversetstring) :: string
        integer(kind=intType) :: i
 
        if (associated(string%nodeData)) &
            deallocate (string%nodeData)
 
        if (associated(string%intNodeData)) &
            deallocate (string%intNodeData)
 
        if (associated(string%conn)) &
            deallocate (string%conn)
 
        if (associated(string%pNodes)) &
            deallocate (string%pNodes)
 
        if (associated(string%pElems)) &
            deallocate (string%pElems)
 
        if (associated(string%cNodes)) &
            deallocate (string%cNodes)
 
        if (associated(string%otherID)) &
            deallocate (string%otherID)
 
        if (associated(string%nte)) &
            deallocate (string%nte)
 
        if (associated(string%subStr)) &
            deallocate (string%subStr)
 
        if (associated(string%elemUsed)) &
            deallocate (string%elemUsed)
 
        if (associated(string%xZipNodeUsed)) &
            deallocate (string%xZipNodeUsed)
 
        if (associated(string%tris)) &
            deallocate (string%tris)
 
        if (associated(string%surfCellID)) &
            deallocate (string%surfCellID)
 
        call nullifystring(string)
 
    end subroutine deallocatestring
 
    subroutine setstringpointers(string)
 
        use constants
        implicit none
        type(oversetstring) :: string
        string%x => string%nodeData(1:3, :)
        string%norm => string%nodeData(4:6, :)
        string%perpNorm => string%nodeData(7:9, :)
        string%h => string%nodeData(10, :)
 
        string%ind => string%intNodeData(1, :)
        string%cluster => string%intNodeData(2, :)
        string%family => string%intNodeData(3, :)
 
    end subroutine setstringpointers
 
    subroutine createorderedstrings(master, strings, nString)
 
        use constants
        implicit none
 
        ! Input/Output
        type(oversetstring) :: master
        type(oversetstring), dimension(:), allocatable :: strings
 
        ! Working
        integer(kind=intType) :: nElems, nNodes, curElem, nString, iStart, i, firstElem
        type(oversetstring), pointer :: stringsLL, str
 
        ! The next step is to create ordered strings based on the
        ! connectivity. This is a purely logical operation. We don't know
        ! how many actual strings we will need so we will use a linked
        ! list as we go.
        call createnodetoelem(master)
 
        ! Allocate some additional arrays we need for doing the chain
        ! searches.
        nelems = master%nElems
        nnodes = master%nNodes
        allocate (master%elemUsed(nelems), master%subStr(2, nelems), &
                  master%cNodes(2, nnodes))
 
        master%cNodes = 0
        master%elemUsed = 0
        curelem = 1
        nstring = 0
        do while (curelem < master%nElems)
 
            ! Arbitrarily get the first node for my element:
            istart = master%conn(1, curelem)
            nelems = master%nte(1, istart)
 
            ! ----------------------
            ! First side of chain:
            ! ----------------------
            firstelem = master%nte(2, istart)
            master%subStr(1, 1) = firstelem
            call dochain(master, istart, 1)
 
            ! ----------------------
            ! Second side of chain:
            ! ----------------------
            if (nelems > 1) then
                firstelem = master%nte(3, istart)
 
                ! Make sure the second one wasn't wrapped around on a
                ! periodic chain
                if (master%elemUsed(firstelem) == 0) then
 
                    master%subStr(2, 1) = firstelem
                    call dochain(master, istart, 2)
                    call combinechainbuffers(master)
                end if
            end if
 
            ! We now have a boundary string stored in master%subString(1,
            ! :nSubStr(1)). These are actually the element numbers of the
            ! master that form a continuous chain.
 
            ! Create or add a new string to our linked list
            ! "stringsLL".
            if (nstring == 0) then
                allocate (stringsll)
                nstring = 1
                stringsll%next => stringsll
                str => stringsll
            else
                allocate (str%next)
                str%next%next => stringsll
                str => str%next
                nstring = nstring + 1
            end if
 
            ! Create a substring from master based on the elements we
            ! have in the buffer
            call createsubstringfromelems(master, str, nstring)
 
            ! Scan through until we find the next unused element:
            do while ((master%elemUsed(curelem) == 1) .and. (curelem < master%nElems))
                curelem = curelem + 1
            end do
        end do
 
        ! Put the strings into an regular array which will be easier to
        ! manipulate.
        allocate (strings(nstring))
        str => stringsll
        i = 0
        do while (i < nstring)
            i = i + 1
            strings(i) = str ! This is derived type assigment.
            call nullifystring(str)
            str => str%next
        end do
 
    end subroutine createorderedstrings
 
    subroutine performselfzip(master, strings, nStrings, debugZipper)
 
        use constants
        use kdtree2_module
        use inputoverset, only: selfzipcutoff
        implicit none
 
        ! Input/Output
        type(oversetstring) :: master
        integer(kind=intType) :: nStrings
        type(oversetstring), dimension(nStrings), target :: strings
        logical, intent(in) :: debugZipper
 
        ! Workging
        type(oversetstring), pointer :: str
        real(kind=realtype) :: cutoff
        integer(kind=intType) :: i, j, nZipped
 
        ! Now determine if there are any "holes" or periodic strings
        ! without anything inside of it. Ie closed loops. If there isn't
        ! we can self zip. Otherwise, we falg it so that it isn't touched
        ! and automatically pocket zipped at th end.
 
        do i = 1, nstrings
            str => strings(i)
            str%isPocket = .true.
            do j = 1, str%nNodes
                if (str%otherID(1, j) /= -1) then
                    str%isPocket = .false.
                end if
            end do
 
            if (.not. str%isPocket) then
                zipperloop: do j = 1, 5
                    if (j == 1) then
                        cutoff = selfzipcutoff
                    else
                        cutoff = 90_realtype
                    end if
                    call selfzip(strings(i), cutoff, nzipped)
                    if (nzipped == 0) then
                        exit zipperloop
                    end if
                end do zipperloop
            end if
        end do
 
        ! Now we need to redo the string matching becuase the self-zip
        ! shortened the strings
        call stringmatch(strings, nstrings, debugzipper)
 
    end subroutine performselfzip
 
    subroutine reducegapstring(string)
 
        ! Generic routine for removing duplicate nodes on the given
        ! string. The string is returned with the nodes and connectivities
        ! adjusted accordingly.
 
        use constants
        use utils, only: pointreduce, mynorm2
        implicit none
 
        ! Input/Ouput
        type(oversetstring), intent(inout) :: string
 
        ! Working:
        real(kind=realtype) :: minedge
        integer(kind=intType) :: nUnqiue, i, n1, n2, nUnique, idx
        integer(kind=intType), dimension(:), allocatable :: link
        real(kind=realtype), dimension(:, :), allocatable :: uniquenodes
        real(kind=realtype), dimension(:, :), pointer :: nodedataptr
        integer(kind=intType), dimension(:, :), pointer :: intNodeDataPtr
        integer(kind=intType), dimension(:), allocatable :: normCounter
        real(kind=realtype), dimension(:, :), allocatable :: uniquenorms
 
        ! We will do a sort of adaptive tolernace here: Get the minium edge
        ! length and base the tolerance on that:
 
        minedge = huge(1.0d0)
 
        do i = 1, string%nElems
            n1 = string%conn(1, i)
            n2 = string%conn(2, i)
            minedge = min(minedge, mynorm2(string%x(:, n1) - string%x(:, n2)))
        end do
 
        allocate (link(string%nNodes), uniquenodes(3, string%nNodes))
 
        call pointreduce(string%x, string%nNodes, minedge / 1000.0, uniquenodes, link, nunique)
 
        ! Now average the normals for any duplicate nodes. This is to handle any discrepancies
        ! for h-type mesh topologies where the surface block is fully represented by the volume connectivity
        allocate (normcounter(nunique), uniquenorms(3, nunique))
        normcounter(:) = zero
        uniquenorms(:, :) = zero
        ! sum the norms for the unique node and count how many duplicates there are for a given node
        do i = 1, string%nNodes
            idx = link(i)
            uniquenorms(:, idx) = uniquenorms(:, idx) + string%nodeData(4:6, i)
            normcounter(idx) = normcounter(idx) + 1
        end do
 
        ! Now divide to get the average and assign back to original data storage
        do i = 1, string%nNodes
            idx = link(i)
            string%nodeData(4:6, i) = uniquenorms(:, idx) / normcounter(idx)
        end do
        deallocate (normcounter, uniquenorms)
        ! Averageing is complete
 
        ! Update the connectivity to use the new set of nodes
        do i = 1, string%nElems
            string%conn(1, i) = link(string%conn(1, i))
            string%conn(2, i) = link(string%conn(2, i))
        end do
 
        ! Reallocate the node based data to the correct size. Set pointers
        ! to original data first.
        nodedataptr => string%nodeData
        intnodedataptr => string%intNodeData
        allocate (string%nodeData(10, nunique), string%intNodeData(3, nunique))
 
        ! Reset the pointers
        call setstringpointers(string)
 
        do i = 1, string%nNodes
            string%nodeData(:, link(i)) = nodedataptr(:, i)
            string%intNodeData(:, link(i)) = intnodedataptr(:, i)
        end do
        string%nNodes = nunique
 
        ! deallocate the pointer data which is actually the original data
        deallocate (nodedataptr, intnodedataptr, link, uniquenodes)
 
    end subroutine reducegapstring
 
    recursive subroutine createnodetoelem(string)
 
        ! Produce the inverse of the connectivity...the nodeToElem
        ! array. Each node should point to 1 element (at a
        ! boundary) or two elements for a normal part of a chain.
 
        use constants
        use utils, only: terminate
        implicit none
 
        ! Input/Output
        type(oversetstring) :: string
 
        ! Working
        integer(kind=intType) :: i, j, ii, jj, n(2), m(2), curelem, ndup
        integer(kind=intType), dimension(string%nElems) :: duplicated
        integer(kind=intType), dimension(:, :), pointer :: tmpconn
        logical :: duplicateelement
 
        allocate (string%nte(3, string%nNodes))
        string%nte = 0
        duplicated = 0
 
        do i = 1, string%nElems
            ! Node numbers we're working with:
            n = string%conn(:, i)
 
            ! For each node check which elements (if any) are already
            ! connected. We need to check them again the node numbers n1 and n2
 
            duplicateelement = .false.
            do jj = 1, 2
 
                do j = 1, string%nte(1, n(jj)) ! Loop over the element numbers already here:
                    curelem = string%nte(j + 1, n(jj))
 
                    m = string%conn(:, curelem)
 
                    if (m(1) == n(1) .and. m(2) == n(2)) then
                        duplicateelement = .true.
                    else if (m(1) == n(2) .and. m(2) == n(1)) then
                        ! Element exists, but it is the wrong order...don't
                        ! know what to do with this, probably an error or
                        ! maybe a corner case I haven't thought of.
                        call terminate("makeBoundaryString", "Inconsistent duplicate edge.")
                    end if
                end do
            end do
 
            if (.not. duplicateelement) then
                do jj = 1, 2
                    string%nte(1, n(jj)) = string%nte(1, n(jj)) + 1
                    ii = string%nte(1, n(jj))
                    string%nte(ii + 1, n(jj)) = i
                end do
            else
                ! Well, we've figured out that this element is actually a
                ! duplicate so we'll make a note of that
                duplicated(i) = 1
            end if
        end do
 
        ! If we have duplicated elements, modify the conn to adjust for this.
        ndup = sum(duplicated)
        if (ndup > 0) then
            tmpconn => string%conn
 
            allocate (string%conn(2, string%nElems - ndup))
 
            j = 0
            do i = 1, string%nElems
                if (duplicated(i) == 0) then
                    j = j + 1
                    string%conn(:, j) = tmpconn(:, i)
 
                end if
            end do
 
            ! Set the new number of elements
            string%nElems = string%nElems - ndup
 
            ! Don't forget to deallocate the tmpConn pointer which is
            ! actually the original conn data.
            deallocate (tmpconn)
 
            ! Destroy nte and call myself again to get the final correct nte
            ! without the duplicates.
            deallocate (string%nte)
            call createnodetoelem(string)
        end if
    end subroutine createnodetoelem
 
    subroutine dochain(master, iStart, iSub)
 
        use constants
        implicit none
        ! Input/OUtput
        type(oversetstring) :: master
        integer(kind=intType), intent(in) :: iStart, iSub
 
        ! Working
        integer(Kind=intType) :: i, j, jj, c, n1, n2, curNode, nextNode
        integer(Kind=intType) :: elem1, elem2, curElem, nextElem
        integer(kind=intType) :: N
 
        ! The number of elements in this substring
        n = 1
 
        curnode = istart
 
        chainloop: do
 
            ! Get the currnet element
            curelem = master%subStr(isub, n)
 
            ! Flag the element as used:
            master%elemUsed(curelem) = 1
 
            ! Get the two nodes for the current element:
            n1 = master%conn(1, curelem)
            n2 = master%conn(2, curelem)
 
            if (n1 == curnode) then
                nextnode = n2
            else
                nextnode = n1
            end if
 
            ! Exit condition 1: Next node was our starting node:
            if (nextnode == istart) then
                exit chainloop
            end if
 
            ! Exit condition 2: The next node has only 1 element, (the one
            ! we're currently on) so that means the the chain is finished
            c = master%nte(1, nextnode)
 
            if (c == 1) then
                exit chainloop
 
            else if (c == 2) then
                ! With c=2 this easy, just extract the two elements
                elem1 = master%nte(2, nextnode)
                elem2 = master%nte(3, nextnode)
 
                if (elem1 == curelem) then
                    nextelem = elem2
                else
                    nextelem = elem1
                end if
            end if
 
            ! Now add the "nextElem" to our chain:
            n = n + 1
            master%subStr(isub, n) = nextelem
 
            ! Flag this elemet as being used
            master%elemUsed(nextelem) = 1
 
            ! Finally set the nextNode back to the current node for the next
            ! iteration
            curnode = nextnode
        end do chainloop
        master%nSubStr(isub) = n
    end subroutine dochain
 
    subroutine createsubstringfromelems(p, s, id)
 
        use constants
        implicit none
 
        ! Input/output
        type(oversetstring), target, intent(in) :: p
        type(oversetstring), intent(out) :: s
        integer(kind=intType), intent(in) :: id
 
        ! Working
        integer(kind=intType) :: i, j, n1, n2, k
        integer(kind=intType), dimension(:), allocatable :: nodeUsed
 
        ! First thing we always have to do with a new string is to nullify
        ! all the poitners
        call nullifystring(s)
 
        ! Firstly we can set the number of elements, since we know precisely
        ! what this is:
 
        s%nElems = p%nSubStr(1)
        s%myID = id
 
        ! Next determine the number of nodes. This is done by flagging the
        ! nodes in the parent that are used by 's'
        allocate (nodeused(p%nNodes))
        nodeused = 0
        k = 0
        do i = 1, s%nElems
            n1 = p%conn(1, p%subStr(1, i))
            n2 = p%conn(2, p%subStr(1, i))
            if (nodeused(n1) == 0) then
                k = k + 1
                nodeused(n1) = k
            end if
 
            if (nodeused(n2) == 0) then
                k = k + 1
                nodeused(n2) = k
            end if
        end do
 
        ! We can now set the number of nodes the substring has
        s%nNodes = k
 
        ! The number of nodes will equal the number of elements iff the
        ! string is period. Otherwise we will have 1 more node than element.
 
        if (s%nNodes == s%nElems) then
            s%isPeriodic = .true.
        end if
 
        ! Allocate and set the node and element parent information
        allocate (s%pElems(s%nElems), s%pNodes(s%nNodes))
 
        do i = 1, s%nElems
            s%pElems(i) = p%subStr(1, i)
        end do
 
        ! Now create the pNodes ("link") array such that pNodes(i) points to
        ! the node index in the parent
        j = 0
        do i = 1, p%nNodes
            if (nodeused(i) /= 0) then
                s%pNodes(nodeused(i)) = i
            end if
        end do
 
        ! Set the parent's cNode to point to my nodes
        do i = 1, s%nNodes
            p%cNodes(:, s%pNodes(i)) = (/s%myID, i/)
        end do
 
        ! Now that we know the mapping between by local nodes-based
        ! quantities and the parent, we can allocate and set all the
        ! node-based quantities.
 
        allocate (s%nodeData(10, s%nNodes), s%intNodeData(3, s%nNodes))
 
        ! Set the string pointers
        call setstringpointers(s)
 
        do i = 1, s%nNodes
            s%nodeData(:, i) = p%nodeData(:, s%pNodes(i))
            s%intNodeData(:, i) = p%intNodeData(:, s%pNodes(i))
        end do
 
        ! We can now create the local conn too, *USING THE LOCAL NODE NUMBERS*
        allocate (s%conn(2, s%nElems))
 
        do i = 1, s%nElems
            s%conn(1, i) = nodeused(p%conn(1, s%pElems(i)))
            s%conn(2, i) = nodeused(p%conn(2, s%pElems(i)))
        end do
 
        ! Set the pointer to my parent.
        s%p => p
 
        deallocate (nodeused)
 
        ! Last thing we can do is create the nodeToElem for the substring.
        call createnodetoelem(s)
    end subroutine createsubstringfromelems
 
    subroutine combinechainbuffers(s)
 
        use constants
        implicit none
        type(oversetstring), intent(inout) :: s
        integer(kind=intType) :: N1, N2
 
        n1 = s%nSubStr(1)
        n2 = s%nSubStr(2)
 
        ! First reverse the direction of string 2 of the nodes we found
        s%subStr(2, 1:n2) = s%subStr(2, n2:1:-1)
 
        ! Now String 1 can be tacked on the end of string2
        s%subStr(2, n2 + 1:n2 + n1) = s%subStr(1, 1:n1)
 
        ! And finally copied back to string1
        s%subStr(1, 1:n1 + n2) = s%subStr(2, 1:n1 + n2)
        s%nSubStr(1) = n1 + n2
 
    end subroutine combinechainbuffers
 
    subroutine selfzip(s, cutOff, nZipped)
 
        use constants
        use kdtree2_module
        use utils, only: mynorm2, cross_prod
        implicit none
 
        ! Input/Output
        type(oversetstring), intent(inout), target :: s
        integer(Kind=intType), intent(out) :: nZipped
        real(kind=realtype), intent(in) :: cutoff
 
        ! Working
        integer(kind=intType) :: i, j, k, N, ii, im1, ip1
        logical :: lastNodeZippered, added
        real(kind=realtype), dimension(3) :: v1, v2, norm
        real(kind=realtype) :: coscutoff, costheta, r2, v1nrm, v2nrm
        integer(Kind=intType), dimension(:), allocatable :: nodeMap
        type(kdtree2_result), dimension(:), allocatable :: results
 
        ! Perform self zipping on the supplied string. The string at this
        ! point should be either peroidic or since sinded --- no multiple
        ! loops should be left. Therefore, we can count on the nodes being
        ! in order.
 
        allocate (results(25))
        allocate (nodemap(s%nNodes))
        nodemap = 1
 
        coscutoff = cos(cutoff * pi / 180)
        nzipped = 0
 
        ! Peroidic string starts at node 1, and uses node 'N' as the previous
        ! node. Single chains start at node 2 and only go to the N-1 node.
        if (s%isPeriodic) then
            im1 = s%nNodes
            ii = 1
            ip1 = 2
            n = s%nNodes
        else
            im1 = 1
            ii = 2
            ip1 = 3
            n = s%nNodes - 1
        end if
 
        do while (ii <= n)
 
            ! Peroidic string at end...loop around
            if (s%isPeriodic .and. ii == n) then
                ip1 = 1
            end if
 
            lastnodezippered = .false.
 
            ! Determine the anlge between the vectors
            v1 = s%x(:, ip1) - s%x(:, ii)
            v2 = s%x(:, im1) - s%x(:, ii)
            v1nrm = mynorm2(v1)
            v2nrm = mynorm2(v2)
            call cross_prod(v2, v1, norm)
            norm = norm / mynorm2(norm)
 
            if (dot_product(norm, s%norm(:, ii)) > zero) then
 
                ! the dot product of the im1 and ip1 nodes have to be close
                if (dot_product(s%norm(:, ip1), s%norm(:, im1)) > 0.80) then
 
                    costheta = dot_product(v1, v2) / (v1nrm * v2nrm)
 
                    if (costheta > coscutoff) then
 
                        call addpotentialtriangle(s, im1, ii, ip1, nodemap, &
                                                  results, added)
 
                        if (added) then
                            nzipped = nzipped + 1
                            lastnodezippered = .true.
                        end if
                    end if
                end if
            end if
 
            if (lastnodezippered) then
                ! Skip the next node...we'll get it on the next pass
                ii = ii + 2
                im1 = ii - 1
                ip1 = ii + 1
            else
                ! Just shuffle along
                ii = ii + 1
                ip1 = ii + 1
                im1 = ii - 1
            end if
        end do
 
        ! Now we will modify our string to remove the elements and nodes
        ! that got knocked off due to self zipping. This way the calling
        ! process still sees the same string, it just gets a little
        ! shorter.
 
        call shortenstring(s, nodemap)
        deallocate (results, nodemap)
 
    end subroutine selfzip
 
    subroutine crosszip(str1, N1, N2, str2, N3, N4, debugZipper, failed)
 
        use constants
        use utils, only: mynorm2, cross_prod
 
        implicit none
 
        type(oversetstring), intent(inout) :: str1, str2
        integer(kind=intType) :: N1, N2, N3, N4
        logical :: debugZipper, failed
        ! Working
        type(oversetstring), pointer :: p
        integer(kind=intType) :: stepsA, stepsB, nStepsA, nStepsB
        integer(kind=intType) :: nTriToAdd, ii, i, j, k, A, B, Ap, Bp
        integer(kind=intType) :: aPrev, bPrev
        real(kind=realtype), dimension(3) :: pta, ptb, ptap, ptbp
        !real(kind=realType), dimension(3) :: ptAPrev, ptBPRev
        real(kind=realtype), dimension(3) :: aoff, boff, apoff, bpoff
        real(kind=realtype), dimension(3) :: norma, normb, normap, normbp
        real(kind=realtype), dimension(3) :: perpa, perpb, perpap, perpbp
        !real(kind=realType), dimension(3) :: normAPrev, normBPrev
        !real(kind=realType), dimension(3) :: perpAPrev, perpBPrev
        real(kind=realtype), dimension(3) :: trinorm1, quadnorm1
        real(kind=realtype), dimension(3) :: trinorm2, quadnorm2
        logical :: aValid, bValid, advanceA, aPreferred, area1, area2
        logical :: advanceB
        logical :: changeA, changeB
        logical :: aValidPrev, bValidPrev, advanceAPrev, advanceBPrev
        real(kind=realtype) :: sum1, sum2, h, dpa, dpb
        !am real(kind=realType), parameter :: cutOff = 0.95*3
        real(kind=realtype), parameter :: cutoff = 0.85 * 3
        ! First determine the the total number of triangles we will add
        ! total. It is equal to the total number of triangles on each
        ! string. This will form the index on the do loop.
        failed = .false.
        ! Str1 goes forward
        if (n2 > n1) then
            nstepsa = n2 - n1
        else if (n2 < n1) then
            nstepsa = n2 + str1%nNodes - n1
        else ! N1 == N2
            nstepsa = str1%nElems
        end if
 
        ! Str2 goes backwards
        if (n3 < n4) then
            nstepsb = n3 + str2%nNodes - n4
        else if (n3 > n4) then
            nstepsb = n3 - n4
        else ! N3 == N4
            nstepsb = str2%nElems
        end if
 
        ! Initialize these out of bounds incase something goes very wrong.
        aprev = -1
        bprev = -1
 
        ! The number of steps we've performed in each edge
        stepsa = 0
        stepsb = 0
 
        ! Initialize the front:
        a = n1
        b = n3
        pta = str1%x(:, a)
        ptb = str2%x(:, b)
 
        norma = str1%norm(:, a)
        normb = str2%norm(:, b)
 
        perpa = str1%perpNorm(:, a)
        perpb = str2%perpNorm(:, b)
 
        ap = nextnode(str1, a, .true.)
        bp = nextnode(str2, b, .false.)
        ptap = str1%x(:, ap)
        ptbp = str2%x(:, bp)
        normap = str1%norm(:, ap)
        normbp = str2%norm(:, bp)
        perpap = str1%perpNorm(:, ap)
        perpbp = str2%perpNorm(:, bp)
 
        ! Cross zip nodes N1 to N2 on str1 to nodes N3 to N4 on str2
        ii = 0
        do while (ii < nstepsa + nstepsb)
 
            avalid = .true.
            bvalid = .true.
            ! ---------------------------------------------------------------
            ! Check 1: Point-in-Triangle test: This test considers the
            ! triangle ABA+ and determines if any of the neighbouring points
            ! on either of the two strings is contained inside the
            ! triangle. If the test is positive, A+ must be rejected. The
            ! same test is repeated for B+.
            ! ---------------------------------------------------------------
 
            if (trioverlap(pta, ptb, ptap, str1, a, ap) .or. &
                trioverlap(pta, ptb, ptap, str2, b, b)) then
                avalid = .false.
            end if
 
            if (trioverlap(pta, ptb, ptbp, str1, a, a) .or. &
                trioverlap(pta, ptb, ptbp, str2, b, bp)) then
                bvalid = .false.
            end if
 
            ! ---------------------------------------------------------------
            ! Check 2: Convex quadrilaterl test: This test considers the
            ! quadrilateral ABB+A+ and determines if it is convex. For
            ! connection to point A+ to be valid, the vector areas of
            ! triangles ABA+ and BB+A+ should have the same size. For
            ! connection to B+ to be valid, the vector areas of trianges ABB+
            ! and AB+A+ should be the same sign.  NOTE THAT THIS TEST DOES NOT
            ! ACTUALLY WORK. IT IS 100% INCORRECT!!! THERE ARE CASES WHERE
            ! THE SIGN OF BOTH AREAS ARE OPPOSITE! IT CANNOT BE SAFELY USED.
            ! ---------------------------------------------------------------
 
            ! area1 = positiveTriArea(ptA, ptB, ptAp, normB)
            ! area2 = positiveTriArea(ptB, ptBp, ptAp, normB)
 
            ! if (area1 .neqv. area2) then
            !    aValid = .False.
            ! end if
 
            ! area1 = positiveTriArea(ptA, ptB, ptBp, normA)
            ! area2 = positiveTriArea(ptAp, ptBp, ptA, normA)
 
            ! if (area1 .neqv. area2) then
            !    bValid = .False.
            ! end if
 
            ! Instead, check if the triangle we're going to add has a
            ! positive or negative vector area
 
            area1 = positivetriarea(pta, ptb, ptap, norma)
            if (area1 .eqv. .false.) then
                avalid = .false.
            end if
 
            area2 = positivetriarea(pta, ptb, ptbp, normb)
            if (area2 .eqv. .false.) then
                bvalid = .false.
            end if
 
            ! ---------------------------------------------------------------
            ! Check 3: Prism volume test: Using the surface normals,
            ! "extrude" a prisim in the direction of each surface normal and
            ! find it's volume. It is is not positive, reject the
            ! triangle. Since we don't have the node off wall, we will have
            ! to make do with the normal vectors and average cell size. We
            ! average the cell size and divide by 1000 to give an approximate
            ! offwall distance. Then we use the norm veectors to offset in
            ! that distance to produce the "off" points.
            ! ---------------------------------------------------------------
 
            ! h = quarter*(str1%h(A) + str1%h(Ap) + str2%h(B) + str2%h(Bp)) / 1000
            ! AOff = ptA + normA * h
            ! BOff = ptB + Bnorm * h
            ! ApOff = ptAp + normAp * h
            ! BpOff = ptBp + normBp * h
 
            ! if (prismVol(A, B, Ap, Aoff, Boff, ApOff) < zero) then
            !    aValid = .False.
            ! end if
 
            ! if (prismVol(A, B, Bp, Aoff, Boff, BpOff) < zero) then
            !    bValid = .False.
            ! end if
 
            ! ---------------------------------------------------------------
            ! Check 4: Interpolation stencil test: This one isn't implemented
            ! ---------------------------------------------------------------
 
            ! ---------------------------------------------------------------
            ! Check 5: Surface normal compatibility test. The surface normal
            ! from the triangle should be pointing (mostly) in the same
            ! direction as the normal of the quad that this triangle shares
            ! and edge with. THIS ALSO DOES NOT WORK! What we have to do
            ! instead, is check the normal tri normal against the node
            ! normals it would be using. This is simplier and is vastly
            ! superior.
            ! ---------------------------------------------------------------
 
            call cross_prod(ptb - pta, ptap - pta, trinorm1)
            trinorm1 = trinorm1 / mynorm2(trinorm1)
 
            ! Compute the sum of the dot product of the nodal norms with the triNorm
            sum1 = dot_product(trinorm1, norma) + dot_product(trinorm1, normb) + &
                   dot_product(trinorm1, normap)
 
            call cross_prod(ptb - pta, ptbp - pta, trinorm2)
            trinorm2 = trinorm2 / mynorm2(trinorm2)
 
            sum2 = dot_product(trinorm2, norma) + dot_product(trinorm2, normb) + &
                   dot_product(trinorm2, normbp)
 
            ! Only use this to help pick one if both are still valid:
            if (avalid .and. bvalid .and. dot_product(trinorm1, trinorm2) < 0.8) then
 
                ! Only use this to help pick one if both are still valid:
 
                if (sum1 < cutoff .and. sum2 > cutoff) then
                    avalid = .false.
 
                else if (sum2 < cutoff .and. sum1 > cutoff) then
                    bvalid = .false.
 
                else if (sum1 < cutoff .and. sum2 < cutoff) then
                    ! Both bad. Take the least bad one
                    if (sum1 > sum2) then
                        bvalid = .false.
                    else
                        avalid = .false.
                    end if
                end if
            end if
 
            ! ---------------------------------------------------------------
            ! Check 6: Front angle test: Try to keep the front as close as
            ! possible to the gap edges.
            ! ---------------------------------------------------------------
 
            ! Triangle ABA+. Original implemnetation
            sum1 = abs(vecangle(pta - ptap, ptb - ptap)) + abs(vecangle(ptbp - ptb, ptap - ptb))
            sum2 = abs(vecangle(pta - ptbp, ptb - ptbp)) + abs(vecangle(ptbp - pta, ptap - pta))
 
            if (sum1 > sum2) then
                apreferred = .true.
            else
                apreferred = .false.
            end if
 
            ! ---------------------------------------------------------------
            ! Check 7: End of string test
            ! ---------------------------------------------------------------
 
            if (a == ap) then
                avalid = .false.
                bvalid = .true.
            end if
 
            if (b == bp) then
                bvalid = .false.
                avalid = .true.
            end if
 
            ! ---------------------------------------------------------------
            ! Decide on the triangle we want to take.
            ! ---------------------------------------------------------------
 
            if (avalid .and. .not. bvalid) then
 
                ! We have no choice but to take A+
                call addtri(a, str1, b, str2, ap, str1)
                advancea = .true.
                advanceb = .false.
                ! Flag the nodes as used
                str1%xZipNOdeUsed(a) = 1
                str1%xZipNOdeUsed(ap) = 1
                str2%xZipNOdeUsed(b) = 1
 
            else if (bvalid .and. .not. avalid) then
 
                ! We have no choice but to take B+
 
                call addtri(a, str1, b, str2, bp, str2)
                str1%xZipNOdeUsed(a) = 1
                str2%xZipNOdeUsed(b) = 1
                str2%xZipNOdeUsed(bp) = 1
 
                advancea = .false.
                advanceb = .true.
 
            else if (avalid .and. bvalid) then
 
                ! We could take either. Use the preferred triangle.
                if (apreferred) then
 
                    call addtri(a, str1, b, str2, ap, str1)
 
                    str1%xZipNOdeUsed(a) = 1
                    str1%xZipNOdeUsed(ap) = 1
                    str2%xZipNOdeUsed(b) = 1
 
                    advancea = .true.
                    advanceb = .false.
 
                else
 
                    call addtri(a, str1, b, str2, bp, str2)
                    str1%xZipNOdeUsed(a) = 1
                    str2%xZipNOdeUsed(b) = 1
                    str2%xZipNOdeUsed(bp) = 1
 
                    advancea = .false.
                    advanceb = .true.
 
                end if
 
            else
 
                ! Things are not looking good...but
 
                if (avalidprev .and. bvalidprev) then
                    ! We might be able to save it! The last triangle we added
                    ! was a choice..both were valid, but we picked one
                    ! because it was preferred. Now we know the one we did
                    ! pick screwed us for the next triangle...go back and
                    ! pick the other one instead!
 
                    ! First 'delete' the triangle by decrementing the tri
                    ! counter and edge counters
                    p => str1%p
                    p%nTris = p%nTris - 1
                    p%nEdges = p%nEdges - 3
 
                    ! Now we determine which one was actually added and add
                    ! the other one instead
 
                    if (advanceaprev) then
                        ! We need to add the old B triangle instead, which
                        ! means the A triangle we had added was bad
                        avalidprev = .false.
                        stepsb = stepsb + 1
                        stepsa = stepsa - 1
 
                        call addtri(aprev, str1, b, str2, bp, str2)
 
                        ! Reset the 'A' data by shuffling backwards: The 'A'
                        ! data is copied to 'Ap' and the 'A' data is restored from Aprev
 
                        ap = a
                        ptap = pta
                        normap = norma
                        perpap = perpa
 
                        a = aprev
                        pta = str1%x(:, a)
                        norma = str1%norm(:, a)
                        perpa = str1%perpNorm(:, a)
 
                        ! Increment the 'B' data since we actually used B
 
                        b = bp
                        ptb = ptbp
                        normb = normbp
                        perpb = perpbp
 
                        ! And get the new data for Bp
                        bp = nextnode(str2, b, .false.)
                        ptbp = str2%x(:, bp)
                        normbp = str2%norm(:, bp)
                        perpbp = str2%perpNorm(:, bp)
 
                        ! We *actually* advanced B so..
                        advancebprev = .true.
                        advanceaprev = .false.
                    else
 
                        ! We need to add the old A triangle, which means the B
                        ! triangle we had added was bad
                        bvalidprev = .false.
                        stepsb = stepsb - 1
                        stepsa = stepsa + 1
                        call addtri(a, str1, bprev, str2, ap, str1)
 
                        ! Reset the 'B' data by shuffling backwards: The 'B'
                        ! data is copied to 'Bp' and the 'B' data is restored from Bprev
 
                        bp = b
                        ptbp = ptb
                        normbp = normb
                        perpbp = perpb
 
                        b = bprev
                        ptb = str2%x(:, b)
                        normb = str2%norm(:, b)
                        perpb = str2%perpNorm(:, b)
 
                        ! Increment the 'A' data since we actually used A
 
                        a = ap
                        pta = ptap
                        norma = normap
                        perpa = perpap
 
                        ! And get the new data for Ap
                        ap = nextnode(str2, a, .true.)
                        ptap = str1%x(:, ap)
                        normap = str1%norm(:, ap)
                        perpap = str1%perpNorm(:, ap)
                        ! We *actually* advanced A so..
                        advanceaprev = .true.
                        advancebprev = .false.
 
                    end if
 
                    ! We *don't* increment ii since this is in essence still
                    ! the "last" iteration. We just cycle and try the current
                    ! one again.
                    if (debugzipper) then
                        print *, 'Saved cross zip from bad front.'
                    end if
                    cycle
 
                end if
 
                advancea = .false.
                advanceb = .false.
 
                ! Ewww. neither triangle is valid. Do not add the triangle
                ! just return and let the cross zip restart. This should
                ! skip over the bad area and the pocket zip can do the bad
                ! region.
                failed = .true.
                return
 
            end if
 
            ! Now we have to shuffle along the string.
            if (advancea .and. .not. advanceb) then
 
                stepsa = stepsa + 1
 
                ! Save a copy of the previous A info
                aprev = a
 
                ! Copy the Ap to A
                a = ap
                pta = ptap
                norma = normap
                perpa = perpap
 
                ! And get the new data for Ap
                ap = nextnode(str1, a, .true.)
                ptap = str1%x(:, ap)
                normap = str1%norm(:, ap)
                perpap = str1%perpNorm(:, ap)
 
            else if (advanceb .and. .not. advancea) then
 
                stepsb = stepsb + 1
 
                ! Save a copy of the previous B info incase we need it
                bprev = b
 
                ! Copy the Bp to B
                b = bp
                ptb = ptbp
                normb = normbp
                perpb = perpbp
 
                ! And get the new data for Bp
                bp = nextnode(str2, b, .false.)
                ptbp = str2%x(:, bp)
                normbp = str2%norm(:, bp)
                perpbp = str2%perpNorm(:, bp)
            end if
 
            ! Save the prevoius valid triangles and what was advanced
            avalidprev = avalid
            bvalidprev = bvalid
            advanceaprev = advancea
            advancebprev = advanceb
 
            ! Finally increment the number of triangles we've used so far.
            ii = ii + 1
        end do
 
    contains
 
        function nextnode(str, i, pos)
 
            implicit none
            type(oversetstring), intent(iN) :: str
            integer(kind=intType), intent(in) :: i
            logical, intent(in) :: pos
            integer(kind=intType) :: nextnode
 
            if (pos) then
                if (stepsa == nstepsa) then
                    nextnode = i
                else
                    nextnode = i + 1
                    if (nextnode > str%nNodes) then
                        if (str%isPeriodic) then
                            ! Loop back around
                            nextnode = 1
                        else
                            ! Leave it at the same node
                            nextnode = i
                        end if
                    end if
                end if
            else
                if (stepsb == nstepsb) then
                    nextnode = i
                else
                    nextnode = i - 1
                    if (nextnode < 1) then
                        if (str%isPeriodic) then
                            ! Loop back around
                            nextnode = str%nNodes
                        else
                            ! Leave it at the same node
                            nextnode = i
                        end if
                    end if
                end if
            end if
        end function nextnode
 
        function vecangle(vec1, vec2)
            implicit none
 
            ! Input/Output
            real(kind=realtype), dimension(3), intent(in) :: vec1, vec2
            real(kind=realtype) :: vecangle
 
            ! Working
            real(kind=realtype), dimension(3) :: veca, vecb
 
            veca = vec1 / mynorm2(vec1)
            vecb = vec2 / mynorm2(vec2)
 
            vecangle = acos(dot_product(veca, vecb))
 
        end function vecangle
 
        function elembetweennodes(str, a, b)
            implicit none
 
            ! Input/Output
            type(oversetstring), intent(in) :: str
            integer(kind=intType), intent(in) :: a, b
            integer(kind=intType) :: elembetweennodes
 
            ! Working
            integer(kind=intType) :: e1, e2, e3, e4
 
            if (str%nte(1, a) == 1) then
                e1 = str%nte(2, a)
                e2 = e1
            else
                e1 = str%nte(2, a)
                e2 = str%nte(3, a)
            end if
 
            if (str%nte(1, b) == 1) then
                e3 = str%nte(2, b)
                e4 = e3
            else
                e3 = str%nte(2, b)
                e4 = str%nte(3, b)
            end if
 
            ! Two of the edges are the same. And this is the one that must
            ! be between the two nodes.
 
            if (e1 == e3 .or. e1 == e4) then
                elembetweennodes = e1
            else
                elembetweennodes = e2
            end if
 
        end function elembetweennodes
 
        function triarea(pt1, pt2, pt3)
 
            use constants
            use utils, only: mynorm2, cross_prod
            implicit none
 
            ! Input/Output
            real(kind=realtype), intent(in), dimension(3) :: pt1, pt2, pt3
            real(kind=realtype) :: triarea
 
            ! Working
            real(kind=realtype), dimension(3) :: norm
 
            call cross_prod(pt2 - pt1, pt3 - pt1, norm)
            triarea = half * mynorm2(norm)
 
        end function triarea
 
    end subroutine crosszip
 
    subroutine addtri(A, sA, B, sB, C, sC)
 
        ! Form a triangle from index 'A' on string 'sA' , index 'B' on
        ! string 'sB' and index 'C' on string 'sC'
 
        use constants
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: a, b, c
        type(oversetstring), intent(in) :: sA, sB, sC
 
        ! Working
        type(oversetstring), pointer :: p
        integer(kind=intType) :: mn1, mn2, mn3
        p => sa%p
 
        p%nTris = p%nTris + 1
 
        ! mn = master node
        mn1 = sa%pNodes(a)
        mn2 = sb%pNodes(b)
        mn3 = sc%pNodes(c)
 
        p%tris(:, p%nTris) = (/mn1, mn2, mn3/)
 
        ! Add these three edges to master list of edges
 
        ! Edge 1:
        p%nEdges = p%nEdges + 1
        p%edges(p%nEdges)%n1 = mn1
        p%edges(p%nEdges)%n2 = mn2
 
        ! Edge 2:
        p%nEdges = p%nEdges + 1
        p%edges(p%nEdges)%n1 = mn2
        p%edges(p%nEdges)%n2 = mn3
 
        ! Edge 3:
        p%nEdges = p%nEdges + 1
        p%edges(p%nEdges)%n1 = mn3
        p%edges(p%nEdges)%n2 = mn1
 
    end subroutine addtri
 
    subroutine makecrosszip(p, strings, nStrings, debugZipper)
 
        use constants
        implicit none
 
        ! Input/output
        integer(kind=intType), intent(in) :: nStrings
        type(oversetstring), intent(inout), target :: p, strings(nStrings)
        type(oversetstring), pointer :: s, s1, s2
        logical, intent(in) :: debugZipper
        ! Working
        integer(kind=intType) :: i, iStart, iEnd, jStart, jEnd, iStart_j, iEnd_j
        integer(kind=intType) :: curOtherID, iString, ii, nextI, curIStart, nIElems_j
        integer(kind=intType) :: nIElemsBeg, nJElemsBeg, nElem1, nElem2
        integer(kind=intType) :: iStart_orig, iEnd_orig, jStart_orig, jEnd_orig
        logical :: fullLoop1, fullLoop2, dummy, failed
        ! The purpose of this routine is to determine the ranges on two
        !  paired strings that are continuously paired and suitable for
        !  performing cross zipping.
 
        ! Allocate arrays to keep track of nodes that have already been
        !  used in cross zipping.
        do i = 1, nstrings
            s => strings(i)
            allocate (s%XzipNodeUsed(s%nNodes))
            s%xZipNodeUsed = 0
        end do
 
        strloop: do istring = 1, nstrings
 
            ! Skip strings that were pocekts
            if (strings(istring)%isPocket) then
                cycle
            end if
 
            ! S1 is the curent '1' string we are working with
            s1 => strings(istring)
 
            ! Find the lowest node number that isn't used:
            curistart = startnode(s1)
            do while (curistart > 0)
 
                if (debugzipper) then
                    print *, '------------------------------------------------'
                    print *, 'Starting string ', s1%myid, 'at index ', curistart
                    print *, '------------------------------------------------'
                end if
 
                istart = curistart
                ! Other ID is the string attached at the current pt.
                curotherid = s1%otherID(1, istart)
 
                if (curotherid == -1) then
                    print *, '*************************************************************************'
                    print *, 'Error during makeCrossZip: Point ', istart, 'does not have a matching point'
                    print *, 'Position: ', s1%x(:, istart)
                    print *, '*************************************************************************'
                    stop
                end if
 
                ! S2 is the current '2' string we are working with
                s2 => strings(curotherid)
                jstart = s1%otherID(2, istart)
 
                ! ---------------- s1 increments -------------
                ! The goal is to increment s1 as far as we can go in the
                ! NEGATIVE direction.
                call tracematch(s1, istart, .false., curotherid, iend, fullloop1)
 
                if (.not. fullloop1) then
                    ! Now set iStart to iEnd. Basically we start right at the
                    ! negative end the chain and traverse in the POSITIVE
                    ! direction.
                    istart = iend
                    call tracematch(s1, istart, .true., curotherid, iend, dummy)
                end if
 
                ! Now, iStart -> iEnd (in the positive order) is the maximum
                ! possible extent that s1 could be connected to s1
                ! over. However, s2 may have something to say about that. We
                ! do the same operation for s2. Note that the orders are reversed.
 
                ! ---------------- s2 increments -------------
                call tracematch(s2, jstart, .true., s1%myID, jend, fullloop2)
 
                ! If the first jnode isnt' actually matched to me, like I am
                ! to him. Therefore skip me, and go to the next one.
                if (jstart == jend .and. .not. fullloop2 .and. &
                    s2%otherID(1, jstart) /= s1%myID) then
                    s1%xZipNodeUsed(curistart) = 1
                    curistart = startnode(s1)
                    cycle
                end if
 
                if (.not. fullloop2) then
                    jstart = jend
                    call tracematch(s2, jstart, .false., s1%myID, jend, dummy)
                end if
 
                if ((istart == iend .and. .not. fullloop1) .or. &
                    (jstart == jend .and. .not. fullloop2)) then
                    ! Can't go anywhere. Flag this node and the next.
 
                    s1%xZipNodeUsed(curistart) = 1
                    curistart = startnode(s1)
                    cycle
                end if
 
                if (debugzipper) then
                    print *, 'Initial Range s1:', istart, iend, fullloop1
                    print *, 'Initial Range s2:', jstart, jend, fullloop2
                end if
 
                ! Save the original start/endn locations
                istart_orig = istart
                iend_orig = iend
                jstart_orig = jstart
                jend_orig = jend
 
                if ((istart == iend .and. fullloop1) .and. &
                    (jstart == jend .and. fullloop2)) then
                    ! s1 fully attached to s2
 
                    call closestsymmetricnode(s1, s2, istart, jstart)
                    iend = istart
                    jend = jstart
 
                else if ((istart == iend .and. fullloop1) .and. .not. fullloop2) then
 
                    ! Project jStart and jEnd onto s1
                    istart = s2%otherID(2, jstart)
                    iend = s2%otherID(2, jend)
 
                else if ((jstart == jend .and. fullloop2) .and. .not. fullloop1) then
 
                    ! Project iStart and iEnd onto s2
                    jstart = s1%otherID(2, istart)
                    jend = s1%otherID(2, iend)
 
                else
 
                    ! part of s1 is attached to part of s2
 
                    nielemsbeg = elemsforrange(s1, istart, iend, .true.)
                    njelemsbeg = elemsforrange(s2, jstart, jend, .false.)
 
                    ! This the "projection" of the 'j' string on the 'i'
                    ! string. Basically this is the range the 'j' string
                    ! wants to "use up" on the i string.
                    istart_j = s2%otherID(2, jstart)
                    iend_j = s2%otherID(2, jend)
 
                    ! These could match up with iStart and iEnd or they could
                    ! not. That's what we need to determine here.
 
                    ! Need to check if iStart_j is "larger" than istart. We
                    ! just increment using nextNode to take care of periodic
                    ! boundaries.
                    if (istart_j /= istart) then
                        ! The starting points are different. Increment iStart
                        ! until we find
                        i = istart
                        do ii = 1, nielemsbeg
                            i = nextnode(s1, i)
                            if (i == istart_j) then
                                istart = istart_j
                                exit
                            end if
                        end do
                    end if
 
                    if (iend_j /= iend) then
                        ! The starting points are different. Decrement jEnd
                        ! until we find
                        i = iend
                        do ii = 1, nielemsbeg
                            i = prevnode(s1, i)
                            if (i == iend_j) then
                                iend = iend_j
                                exit
                            end if
                        end do
                    end if
 
                    ! Now with the updated range. Project the iRange back to the
                    ! the final J range.
                    jstart = s1%otherID(2, istart)
                    jend = s1%otherID(2, iend)
 
                end if
 
                if (debugzipper) then
                    print *, 'Zipping string: ', s1%myid, ' with ', s2%myid
                    print *, 's1 range:', istart, iend
                    print *, 's2 range:', jstart, jend
                end if
 
                ! Before we do the zip, make sure the ranges have not
                ! degenerated to 0 elements
                nelem1 = 1 ! Doesn't matter, just not zero
                nelem2 = 1 ! Doesn't matter, just not zero
                if (.not. fullloop1) then
                    nelem1 = elemsforrange(s1, istart, iend, .true.)
                    if (istart == iend) then
                        nelem1 = 0
                    end if
 
                end if
                if (.not. fullloop2) then
                    nelem2 = elemsforrange(s2, jstart, jend, .false.)
                    if (jstart == jend) then
                        nelem2 = 0
                    end if
                end if
 
                if (nelem1 > 0 .and. nelem2 > 0) then
                    ! Do actual cross zip if we still have elements left on both strings
                    call crosszip(s1, istart, iend, s2, jstart, jend, debugzipper, failed)
 
                    ! If we succefully cross zippered what we were suppoed to
                    ! flag all the nodes from the original region as done.
                    if (.not. failed) then
                        call flagnodesused(s1, istart_orig, iend_orig, .true.)
                        call flagnodesused(s2, jstart_orig, jend_orig, .false.)
                    else
                        ! UhOh. We got stopped part way through. Flag just the
                        ! nodes at the beginning that we didn't use. Leave
                        ! those for the pocket. The nodes up to where theh
                        ! cross zip stopped were flagged internally in cross zip.
                        call flagnodesused(s1, istart_orig, istart, .true.)
                        call flagnodesused(s2, jstart_orig, jstart, .false.)
                    end if
                else
 
                    ! Flag the full range of elements are consumed even
                    ! though we didn't do the cross zip. Leave it for the
                    ! pocket zipping.
                    call flagnodesused(s1, istart_orig, iend_orig, .true.)
                    call flagnodesused(s2, jstart_orig, jend_orig, .false.)
                end if
 
                ! Find the next starting index:
                curistart = startnode(s1)
 
            end do
        end do strloop
 
    contains
 
        function startnode(s)
            ! Determine the lowest index of a non-used xzip node for
            ! string 's'.
            implicit none
            type(oversetstring) :: s
            integer(kind=intType) :: startnode, i
 
            ! This will be the return value if all nodes are used:
            startnode = 0
            nodeloop: do i = 1, s%nNodes
                if (s%xZipNodeUsed(i) == 0) then
                    startnode = i
                    exit nodeloop
                end if
            end do nodeloop
        end function startnode
 
        function nextnode(s, i)
 
            implicit none
            type(oversetstring), intent(iN) :: s
            integer(kind=intType), intent(in) :: i
            integer(kind=intType) :: nextnode
 
            ! Normally just increment:
            nextnode = i + 1
 
            if (i == s%nNodes) then
                if (s%isPeriodic) then
                    nextnode = 1
                else
                    ! Can't go any further
                    nextnode = i
                end if
            end if
 
            ! If the next node is used. The next node is set the current
            ! one.
            if (s%xZipNodeUsed(nextnode) == 1) then
                nextnode = i
            end if
        end function nextnode
 
        function simplenextnode(s, i)
 
            implicit none
            type(oversetstring), intent(iN) :: s
            integer(kind=intType), intent(in) :: i
            integer(kind=intType) :: simplenextnode
 
            ! Normally just increment:
            simplenextnode = i + 1
 
            if (i == s%nNodes) then
                if (s%isPeriodic) then
                    simplenextnode = 1
                else
                    ! Can't go any further
                    simplenextnode = i
                end if
            end if
        end function simplenextnode
 
        function prevnode(s, i)
 
            implicit none
            type(oversetstring), intent(iN) :: s
            integer(kind=intType), intent(in) :: i
            integer(kind=intTYpe) :: prevnode
            ! Normally just increment:
            prevnode = i - 1
 
            if (i == 1) then
                if (s%isPeriodic) then
                    prevnode = s%nNodes
                else
                    ! Can't go any further
                    prevnode = i
                end if
            end if
 
            ! If the next node is used. The next node is set the current
            ! one.
            if (s%xZipNodeUsed(prevnode) == 1) then
                prevnode = i
            end if
        end function prevnode
 
        subroutine tracematch(s, iStart, pos, checkID, iEnd, fullLoop)
 
            implicit none
 
            ! Given a starting position 'iStart' on string 's', traverse in
            ! the 'POSitive' or '.not. POSitive' direction checking that the
            ! otherID still matches "checkID". Return the ending position
            ! 'iEnd'.
 
            ! Input/Output
            type(oversetstring) :: s
            integer(kind=intType), intent(in) :: iStart, checkID
            logical, intent(in) :: pos
            integer(kind=intType), intent(out) :: iEnd
            logical, intent(out) :: fullLoop
 
            ! Working
            integer(kind=intType) :: i, nextI
 
            i = istart
            fullloop = .false.
 
            traverseloop: do
                if (pos) then
                    nexti = nextnode(s, i)
                else
                    nexti = prevnode(s, i)
                end if
                if (nexti == i .or. s%otherID(1, nexti) /= checkid) then
                    ! We can't go any further than we already are
                    iend = i
                    exit traverseloop
                end if
 
                ! Continue to the next one.
                i = nexti
 
                if (i == istart) then
                    fullloop = .true.
                    iend = i
                    exit traverseloop
                end if
            end do traverseloop
        end subroutine tracematch
 
        subroutine flagnodesused(s, N1, N2, pos)
 
            implicit none
 
            ! Input/Output
            type(oversetstring) :: s
            integer(kind=intType), intent(in) :: N1, N2
            logical, intent(in) :: pos
 
            ! Working
            integer(kind=intType) :: nSteps, i, nextI
 
            if (pos) then
                if (n2 > n1) then
                    nsteps = n2 - n1
                else if (n2 < n1) then
                    nsteps = n2 + s%nNodes - n1
                else ! N1 == N2
                    nsteps = s%nElems
                end if
            else
                if (n1 < n2) then
                    nsteps = n1 + s%nNodes - n2
                else if (n1 > n2) then
                    nsteps = n1 - n2
                else ! N3 == N4
                    nsteps = s%nElems
                end if
            end if
 
            s%xZipNodeUsed(n1) = 1
            i = n1
            do ii = 1, nsteps
                if (pos) then
                    nexti = nextnode(s, i)
                else
                    nexti = prevnode(s, i)
                end if
 
                s%xZipNodeUsed(nexti) = 1
                i = nexti
            end do
        end subroutine flagnodesused
 
        function elemsforrange(s, N1, N2, pos)
            ! Determine the number of elements between N1 and N2 for for the
            ! "POSitive" or "not POSIitive" (negative) direction.
 
            implicit none
            type(oversetstring) :: s
            integer(kind=intType), intent(in) :: n1, n2
            logical :: pos
            integer(kind=intType) :: elemsforrange
 
            if (.not. s%isPeriodic) then
                if (pos) then
                    elemsforrange = n2 - n1
                else
                    elemsforrange = n1 - n2
                end if
            else ! Periodic
                if (pos) then
                    if (n2 == n1) then
                        elemsforrange = s%nElems
                    else if (n2 > n1) then
                        elemsforrange = n2 - n1
                    else
                        elemsforrange = n2 + s%nNodes - n1
                    end if
                else
                    if (n1 == n2) then
                        elemsforrange = s%nElems
                    else if (n1 > n2) then
                        elemsforrange = n1 - n2
                    else
                        elemsforrange = n1 + s%nNodes - n2
                    end if
                end if
            end if
        end function elemsforrange
 
    end subroutine makecrosszip
 
    subroutine makepocketzip(p, strings, nStrings, pocketMaster, debugZipper)
        use constants
        use oversetdata, only: oversetstring, oversetedge
        use oversetutilities, only: qsortedgetype
        use kdtree2_module
        implicit none
 
        ! Input/output
        integer(kind=intType), intent(in) :: nStrings
        type(oversetstring), intent(in) :: p, strings(nStrings)
        type(oversetstring) :: pocketMaster
        logical, intent(in) :: debugZipper
 
        ! Local variables
        integer(kind=intType) :: i, j, nsum1, nsum2, ndiff1, ndiff2, ipedge, icur
        integer(kind=intType) :: n1, n2, npolyEdges
        integer(kind=intType) :: nNodes1, nNodes2, cn1, cn2, str1, str2
        type(oversetedge), allocatable, dimension(:) :: polyEdges
        type(oversetedge) :: e1, e2
        type(oversetstring), pointer :: stringsLL, str
        integer(kind=intType) :: npocketEdges, nFullStrings, nNodes
        integer(kind=intType) :: ip, curElem, nElems, iStart, firstElem
        type(oversetstring), allocatable, dimension(:), target :: pocketStringsArr
 
        ! ---------------------------------------------------------------
        ! PocketZip 1:
        ! First sort the edges.
        ! ---------------------------------------------------------------
        call qsortedgetype(p%Edges, p%nEdges)
 
        ! Now gather up the left-over edges for pocket zipping.
 
        ! Over estimate of remaining pocket edges to zip
        allocate (polyedges(p%nEdges))
 
        ! Eliminate the edges going through the ordered edges.
        ! The sorted opposite edges are canceled in pairs.
        npolyedges = 0
        i = 1
        do while (i <= p%nEdges)
 
            if (i == p%nEdges) then
                ! This must be the last free edge:
                e1 = p%Edges(i)
                npolyedges = npolyedges + 1
                polyedges(npolyedges) = e1
                i = i + 1
                cycle
            end if
 
            ! Two edges in sequence
            e1 = p%Edges(i)
            e2 = p%Edges(i + 1)
 
            ! First determine if e1 is at the end of two single ended
            ! chains. In this case the edge *will* not be paired and that's
            ! correct.
 
            str1 = p%cNodes(1, e1%n1) ! node1's child fullStrings ID
            cn1 = p%cNodes(2, e1%n1) ! node1's child fullStrings node index
            nnodes1 = strings(str1)%nNodes ! node1's child fullStrings nNodes size
 
            str2 = p%cNodes(1, e1%n2) ! node2's child fullStrings ID
            cn2 = p%cNodes(2, e1%n2) ! node2's child fullStrings node index
            nnodes2 = strings(str2)%nNodes ! node2's child fullStrings nNodes size
 
            if (str1 /= str2) then
                if (.not. strings(str1)%isperiodic .and. &
                    .not. strings(str2)%isperiodic .and. &
                    (cn1 == 1 .or. cn1 == nnodes1) .and. (cn2 == 1 .or. cn2 == nnodes2)) then
                    ! Increment just 1 in 1 to skip over edge e1.
                    i = i + 1
                    cycle
                end if
            end if
 
            ! The sum and difference:
            nsum1 = e1%n1 + e1%n2
            nsum2 = e2%n1 + e2%n2
 
            ndiff1 = e1%n2 - e1%n1
            ndiff2 = e2%n2 - e2%n1
 
            if (nsum1 == nsum2 .and. ndiff1 + ndiff2 == 0) then
                ! These edges cancel. Great.
                i = i + 2
                cycle
            else
                ! Add just the first edge
                npolyedges = npolyedges + 1
                polyedges(npolyedges) = e1
                i = i + 1
            end if
        end do
 
        ! Define pocketMaster string
        call nullifystring(pocketmaster)
        pocketmaster%myID = 88
        pocketmaster%nElems = npolyedges
        pocketmaster%nNodes = npolyedges * 2
        pocketmaster%nEdges = 0
        allocate (pocketmaster%nodeData(10, 2 * npolyedges), &
                  pocketmaster%intNodeData(3, 2 * npolyedges), &
                  pocketmaster%conn(2, npolyedges))
 
        ! Dump the data into the pocketMaster
        do i = 1, npolyedges
            pocketmaster%nodeData(:, 2 * i - 1) = p%nodeData(:, polyedges(i)%n1)
            pocketmaster%intNodeData(:, 2 * i - 1) = p%intNodeData(:, polyedges(i)%n1)
 
            pocketmaster%nodeData(:, 2 * i) = p%nodeData(:, polyedges(i)%n2)
            pocketmaster%intNodeData(:, 2 * i) = p%intNodeData(:, polyedges(i)%n2)
            pocketmaster%conn(:, i) = (/2 * i, 2 * i - 1/)
        end do
 
        call setstringpointers(pocketmaster)
        call reducegapstring(pocketmaster)
        call createnodetoelem(pocketmaster)
 
        ! The next step is to create ordered strings based on the
        ! connectivity. This is a purely logical operation. We don't know
        ! how many actual strings we will need so we will use a linked
        ! list as we go.
 
        ! Allocate some additional arrays we need for doing the chain
        ! searches.
        nelems = pocketmaster%nElems
        nnodes = pocketmaster%nNodes
        allocate (pocketmaster%elemUsed(nelems), pocketmaster%subStr(2, nelems), &
                  pocketmaster%cNodes(2, nnodes))
 
        pocketmaster%cNodes = 0
        pocketmaster%elemUsed = 0
        curelem = 1
        nfullstrings = 0
 
        do while (curelem < pocketmaster%nElems)
 
            ! Arbitrarily get the first node for my element:
            istart = pocketmaster%conn(1, curelem)
            nelems = pocketmaster%nte(1, istart)
 
            firstelem = pocketmaster%nte(2, istart)
            pocketmaster%subStr(1, 1) = firstelem
            call dochain(pocketmaster, istart, 1)
 
            ! We now have a boundary string stored in master%subString(1,
            ! :nSubStr(1)). These are actually the element numbers of the
            ! master that form a continuous chain.
 
            ! Create or add a new string to our linked list
            ! "stringsLL".
            if (nfullstrings == 0) then
                allocate (stringsll)
                nfullstrings = 1
                stringsll%next => stringsll
                str => stringsll
            else
                allocate (str%next)
                str%next%next => stringsll
                str => str%next
                nfullstrings = nfullstrings + 1
            end if
 
            ! Create a substring from master based on the elements we
            ! have in the buffer
            call createsubstringfromelems(pocketmaster, str, nfullstrings)
 
            ! Scan through until we find the next unused element:
            do while ((pocketmaster%elemUsed(curelem) == 1) .and. &
                      (curelem < pocketmaster%nElems))
                curelem = curelem + 1
            end do
        end do
 
        ! Temporary strings array for plotting and pocketZipping
        allocate (pocketstringsarr(nfullstrings))
        str => stringsll
        i = 0
        do while (i < nfullstrings)
            i = i + 1
            pocketstringsarr(i) = str ! Derived type assignment
            call nullifystring(str)
            str => str%next
        end do
 
        ! Allocate space for pocket triangles.
        ! (n-sided polygon -> n-2 triangles)
        allocate (pocketmaster%tris(3, 10 * pocketmaster%nElems))
        allocate (pocketmaster%edges(4 * pocketmaster%nElems))
        pocketmaster%nTris = 0
 
        ! Build the pocketMaster tree
        pocketmaster%tree => kdtree2_create(pocketmaster%x, sort=.true.)
 
        if (debugzipper) then
            open (unit=101, file="strings_pocket.dat", form='formatted')
            write (101, *) 'TITLE = "PocketStrings Data" '
 
            write (101, *) 'Variables = "X" "Y" "Z" "Nx" "Ny" "Nz" "Vx" "Vy" "Vz" "ind" &
                 &"gapID" "gapIndex" "otherID" "otherIndex" "ratio"'
            do i = 1, nfullstrings
                ! Temporarily allocate otherID
                allocate (pocketstringsarr(i)%otherID(2, pocketstringsarr(i)%nNodes))
                pocketstringsarr(i)%otherID = -1
 
                call writeoversetstring(pocketstringsarr(i), pocketstringsarr, &
                                        nfullstrings, 101)
            end do
            close (101)
        end if
 
        ! Loop over pocketStrings and begin pocketZip starting
        ! from smallest convex ear.
        do i = 1, nfullstrings
            if (debugzipper) then
                print *, 'Pocket Zipping String ', i, ' of ', nfullstrings
            end if
            pocketziploop: do while (pocketstringsarr(i)%nNodes > 2)
                ! Each pass zips one triangle. Keep zipping
                ! until last triangle is zipped in the pocket polygon.
                call pocketzip(pocketstringsarr(i))
            end do pocketziploop
        end do
 
        ! Destroy the strings array
        do i = 1, nfullstrings
            call deallocatestring(pocketstringsarr(i))
        end do
        deallocate (pocketstringsarr, polyedges)
 
    end subroutine makepocketzip
 
    subroutine pocketzip(s)
 
        use constants
        use kdtree2_module
        use utils, only: mynorm2, cross_prod
 
        implicit none
 
        ! Input parameters
        type(oversetstring), intent(inout), target :: s
 
        ! Local variables
        integer(kind=intType) :: i, j, k, ii, im1, ip1, N
        integer(kind=intType) :: nNodes, nElems, iimin
        real(kind=realtype), dimension(3) :: v1, v2, norm, c
        real(kind=realtype) :: coscutoff, costheta, r2, v1nrm, v2nrm, costhetamax
        real(kind=realtype) :: dp, dpmax
 
        integer(Kind=intType), dimension(:), allocatable :: nodeMap, badNode
        type(kdtree2_result), dimension(:), allocatable :: results
        logical :: added, iiMinSet
        real(kind=realtype), parameter :: fact = 0.95_realtype
        n = s%nNodes
        allocate (results(25), nodemap(n), badnode(n))
        nodemap = 1
        badnode = 0 ! Will become 1 if bad
        outerziploop: do
 
            ! No choice for the last triangle:
            if (n == 3) then
                ii = 1
                im1 = prevnode(ii)
                ip1 = nextnode(ii)
                ! We don't call addPotentialTriangle because we don't have a
                ! choice anymore. Just call the raw addTri command
                call addtri(ip1, s, ii, s, im1, s)
                ! and flag the node as gone
                nodemap(ii) = 0
                exit outerziploop
            end if
 
            iiminset = .false.
 
            ! First find the largest dot product:
            dpmax = -one
            nodeloop1: do ii = 1, n
 
                if (badnode(ii) == 1) then
                    cycle nodeloop1
                end if
 
                ip1 = nextnode(ii)
                im1 = prevnode(ii)
 
                ! Determine the angle between the vectors
                v1 = s%x(:, im1) - s%x(:, ii)
                v2 = s%x(:, ip1) - s%x(:, ii)
                v1nrm = mynorm2(v1)
                v2nrm = mynorm2(v2)
                call cross_prod(v2, v1, norm)
                norm = norm / mynorm2(norm)
                dpmax = max(dpmax, dot_product(norm, s%norm(:, ii)))
            end do nodeloop1
 
            ! Next find the largest cosTheta that is winthin a factor
            !  of dpMax
            costhetamax = -large
            nodeloop2: do ii = 1, n
 
                if (badnode(ii) == 1) then
                    cycle nodeloop2
                end if
 
                ip1 = nextnode(ii)
                im1 = prevnode(ii)
 
                ! Determine the angle between the vectors
                v1 = s%x(:, im1) - s%x(:, ii)
                v2 = s%x(:, ip1) - s%x(:, ii)
                v1nrm = mynorm2(v1)
                v2nrm = mynorm2(v2)
                call cross_prod(v2, v1, norm)
                norm = norm / mynorm2(norm)
                dp = dot_product(norm, s%norm(:, ii))
                if (dp > dpmax * fact) then ! We take this
                    costheta = dot_product(v1, v2) / (v1nrm * v2nrm)
                    if (costheta > costhetamax) then
                        costhetamax = costheta
                        iiminset = .true.
                        iimin = ii
                    end if
                end if
            end do nodeloop2
 
            if (iiminset) then
                ! Zip about node "iimin" if it was set:
                ii = iimin
                ip1 = nextnode(ii)
                im1 = prevnode(ii)
                call addpotentialtriangle(s, ip1, ii, im1, nodemap, results, added)
                if (added) then
                    ! This triangle was good!
                    exit outerziploop
                else
                    ! Bad node. Need to cycle through rest of pocket nodes.
                    ! Remember this bad node in next cycle.
                    badnode(ii) = 1
                    cycle outerziploop
                end if
            else
                ! What does this mean? We didn't find any node to zip. Are they all bad?
                print *, 'Problem with pocket zipper. Somehow we were not able to find "&
                     &"node to add a triangle on. This should not happen. Contact the "&
                     &"Developers.'
                stop
            end if
        end do outerziploop
 
        ! Modify the pocketStrings to remove the two elements and the node
        ! that got eliminated due to pocketZipping.
        call shortenstring(s, nodemap)
        deallocate (nodemap, badnode, results)
 
    contains
        function nextnode(ii)
            implicit none
            integer(kind=intType) :: ii, nextnode
            nextnode = ii + 1
            if (ii == n) then
                nextnode = 1
            end if
        end function nextnode
 
        function prevnode(ii)
            implicit none
            integer(kind=intType) :: ii, prevnode
            prevnode = ii - 1
            if (ii == 1) then
                prevnode = n
            end if
        end function prevnode
    end subroutine pocketzip
 
    subroutine computetrisurfarea(master, area)
 
        ! Computes area sum of all triangles belonging to object master
        use constants
        use utils, only: mynorm2, cross_prod
        implicit none
 
        ! Input parameters
        type(oversetstring), intent(in) :: master
        real(kind=realtype), intent(out) :: area
 
        ! Local variables
        integer(kind=intType) :: i, n1, n2, n3
        real(kind=realtype), dimension(3) :: v1, v2, norm
 
        area = 0.0
        do i = 1, master%nTris
            n1 = master%tris(1, i)
            n2 = master%tris(2, i)
            n3 = master%tris(3, i)
 
            v1 = master%x(:, n2) - master%x(:, n1)
            v2 = master%x(:, n3) - master%x(:, n1)
            call cross_prod(v1, v2, norm)
            area = area + half * mynorm2(norm)
        end do
 
    end subroutine computetrisurfarea
 
    function trioverlap(pt1, pt2, pt3, str, i1, i2)
 
        use constants
        use utils, only: mynorm2, cross_prod
        implicit none
 
        ! Input/Output
        real(kind=realtype), dimension(3), intent(in) :: pt1, pt2, pt3
        integer(kind=intType), intent(in) :: i1, i2
        type(oversetstring), intent(in) :: str
 
        ! Working
        logical :: trioverlap, intri
        integer(kind=intType) :: i
        real(kind=realtype) :: trinorm(3)
 
        ! Note: This is a dumb loop. We need to do a spatial serch here to
        ! only check the nodes around the current point.
 
        call cross_prod(pt2 - pt1, pt3 - pt1, trinorm)
        trinorm = trinorm / mynorm2(trinorm)
 
        trioverlap = .false.
        do i = 1, str%nNodes
            if (i /= i1 .and. i /= i2) then
                if (dot_product(str%norm(:, i), trinorm) > 0.8) then
                    call pointintriangle(pt1, pt2, pt3, str%x(:, i), intri)
                    if (intri) then
                        trioverlap = .true.
                        exit
                    end if
                end if
            end if
        end do
    end function trioverlap
 
    subroutine shortenstring(s, nodeMap)
 
        ! This is an auxilary routine that take a string 's', and a node
        ! map of len s%nNodes, with 1 or 0. A 1 means that the node will
        ! be in the shortened string, 0 means that the node should be
        ! deleted.
        use constants
        implicit none
 
        ! Input/Output
        type(oversetstring) :: s
        integer(kind=intType), dimension(:), intent(inout) :: nodemap
 
        ! Working
        integer(kind=intType) :: nnodes, nelems, nremoved, i, j
        real(kind=realtype), dimension(:, :), pointer :: nodedatatmp
        integer(kind=intType), dimension(:, :), pointer :: conntmp, intnodedatatmp
        integer(kind=intType), dimension(:), pointer :: pnodestmp
 
        ! Now we will modify our string to remove the elements and nodes
        ! that got knocked off due to self zipping. This way the calling
        ! process still sees the same string, it just gets a little
        ! shorter.
 
        ! Save pointers to existing data
        nnodes = s%nNodes
        nelems = s%nElems
        nodedatatmp => s%nodeData
        intnodedatatmp => s%intNodeData
        conntmp => s%conn
        pnodestmp => s%pNodes
 
        ! Convert the nodeMap which currently contains a one if the node
        ! still exists and 0 if it doesn't. This will convert it to the new
        ! node numbers. Ie nodeMap(i) gives the new node index of the
        ! shorted chain. If nodeMap(i) = 0, it is no longer part of the
        ! chain.
        j = 0
        nremoved = 0
        do i = 1, s%nNodes
            if (nodemap(i) == 1) then
                j = j + 1
                nodemap(i) = j
            else
                nremoved = nremoved + 1
            end if
        end do
 
        !  Update the cNodes in the parent so they point to the updated node
        ! numbers. Note that the nodes that have been eliminated, have cNode
        ! = 0, which will identify that it no longer has a child node.
        do i = 1, s%nNodes
            s%p%cNodes(:, s%pNodes(i)) = (/s%myID, nodemap(i)/)
        end do
 
        ! Update the number of nodes/elems in our shorted chain. Every
        ! zipper reduces the number of nodes and number of elems by 1
        s%nNodes = s%nNodes - nremoved
        s%nElems = s%nElems - nremoved
 
        allocate (s%nodeData(10, s%nNodes), s%intNodeData(3, s%nNodes), &
                  s%pNodes(s%nNodes), s%conn(2, s%nElems))
 
        ! Set the pointers for the new string
        call setstringpointers(s)
 
        do i = 1, nnodes
            if (nodemap(i) /= 0) then
                s%nodeData(:, nodemap(i)) = nodedatatmp(:, i)
                s%intNodeData(:, nodemap(i)) = intnodedatatmp(:, i)
                s%pNodes(nodemap(i)) = pnodestmp(i)
            end if
        end do
 
        ! Since we know the string was in order, we can simply redo the connectivity
        do i = 1, s%nElems
            s%conn(:, i) = (/i, i + 1/)
        end do
 
        if (s%isPeriodic) then
            s%conn(2, s%nElems) = 1
        end if
 
        ! Dellocate the existing memory
        deallocate (nodedatatmp, intnodedatatmp, conntmp, pnodestmp)
 
        ! Recreate the node to elem
        if (s%nNodes >= 3) then
            call createnodetoelem(s)
        end if
 
    end subroutine shortenstring
 
    subroutine addpotentialtriangle(s, im1, ii, ip1, nodeMap, results, added)
 
        ! Common routine (for pocketZip and selfZip) to potentially add a
        ! triangle resulting from a single string.
        use constants
        use kdtree2_priority_queue_module
        use kdtree2_module
        implicit none
 
        ! Input/Output
        type(oversetstring) :: s
        integer(kind=intType), intent(in) :: im1, ii, ip1
        integer(kind=intType), intent(inout), dimension(:) :: nodeMap
        type(kdtree2_result), dimension(:), allocatable :: results
        logical, intent(out) :: added
        ! Working:
        real(kind=realtype) :: r2
        real(kind=realtype), dimension(3) :: v1, v2, norm, c
        integer(kind=intType) :: nFound, nalloc, idx, k, j, i
        logical :: overlapFound, inTri
 
        ! We may have a valid triangle. We need to make sure we
        ! don't overlap anyone else.
        !
        ! xim1 +
        !      | \
        !      |   \
        !      |     c
        !      |       \
        !      |         \
        !      +----------+
        !      xi         xip1
        ! We do a ball search based at 'c' which is just the
        ! (average of xip1 and xim1) using a radius defined as the
        ! maximum of (the distance between 'c' and 'xi', half
        ! length of xip1 to xim1)
        !
        added = .false.
 
        c = half * (s%x(:, ip1) + s%x(:, im1))
        r2 = (c(1) - s%x(1, ii))**2 + (c(2) - s%x(2, ii))**2 + (c(3) - s%x(3, ii))**2
 
        r2 = max(r2, (s%x(1, ip1) - s%x(1, im1))**2 + (s%x(2, ip1) - s%x(2, im1))**2 + &
                 (s%x(3, ip1) - s%x(3, im1))**2)
 
        nfound = 0
        outerloop: do
            nalloc = size(results)
            call kdtree2_r_nearest(s%p%tree, c, r2, nfound, nalloc, results)
            if (nfound < nalloc) then
                exit outerloop
            end if
 
            ! Allocate more space and keep going
            deallocate (results)
            nalloc = nalloc * 2
            allocate (results(nalloc))
        end do outerloop
 
        ! We can now be sure that we have all the points inside our
        ! ball. Next we proceed to systematically check them.
        overlapfound = .false.
        nodefoundloop: do k = 1, nfound
            ! Note that we do check nodes from our own string,
            ! except for the the three nodes we're dealing
            ! with. Remember that we are working in our parent's
            ! ording here.
            idx = results(k)%idx
 
            notpartoftriangle: if (idx /= s%pNodes(im1) .and. &
                                   idx /= s%pNodes(ii) .and. idx /= s%pNodes(ip1)) then
 
                ! Only check if the node normal of the point we're
                ! checking is in the same direction as the triangle.
                if (dot_product(s%norm(:, ii), s%p%norm(:, idx)) > zero) then
 
                    ! Finally do the actual trianlge test
                    call pointintriangle(s%x(:, ip1), s%x(:, ii), s%x(:, im1), &
                                         s%p%x(:, idx), intri)
                    if (intri) then
                        ! As soon as 1 is in the triangle, we know the
                        ! triangle is no good.
                        overlapfound = .true.
                        exit nodefoundloop
                    end if
                end if
            end if notpartoftriangle
        end do nodefoundloop
 
        if (.not. overlapfound) then
 
            ! This triangle is good!
            added = .true.
 
            ! Call the generic addTri Routine. Here all the ndoes from the
            ! triangle come from the same string.
            call addtri(ip1, s, ii, s, im1, s)
 
            ! Flag this node as gone
            nodemap(ii) = 0
 
        end if
    end subroutine addpotentialtriangle
 
    subroutine closestsymmetricnode(s1, s2, i, j)
        use constants
        use utils, only: mynorm2
 
        implicit none
 
        ! Input/Output
        type(oversetstring) :: s1, s2
        integer(kind=intType), intent(out) :: i, j
        real(kind=realtype) :: mindist, dist
        integer(kind=intType) :: ii
 
        ! Working:
        mindist = large
 
        do ii = 1, s1%nNodes
            ! "The other index of the matching node on the other string is
            ! me" ie. "I point to you and you point to me"
 
            if (s2%otherID(2, s1%otherID(2, ii)) == ii) then
 
                dist = mynorm2(s1%x(:, ii) - s2%x(:, s1%otherID(2, ii)))
 
                if (dist < mindist) then
                    mindist = dist
                    i = ii
                    j = s1%otherID(2, ii)
                end if
            end if
        end do
 
    end subroutine closestsymmetricnode
 
    subroutine stringmatch(strings, nStrings, debugZipper)
        use constants
        use kdtree2_priority_queue_module
        use kdtree2_module
        use utils, only: mynorm2
 
        implicit none
 
        ! Input/output
        integer(kind=intType), intent(in) :: nStrings
        type(oversetstring), dimension(nstrings), target :: strings
        logical, intent(in) :: debugZipper
 
        ! Working
        integer(kind=intType) :: i, j, k, idx, oid(4)
        integer(kind=intType) :: nAlloc, nUnique, nSearch
        type(kdtree2_result), allocatable, dimension(:) :: results
        type(oversetstring), pointer :: str, master
        logical :: checkLeft, checkRight, concave
        logical :: checkLeft2, checkRight2, concave2
        logical :: leftOK, rightOK, isEndNode
        real(kind=realtype), dimension(3) :: xj, xjp1, xjm1, normj
        real(kind=realtype), dimension(3) :: xk, xkp1, xkm1, normk
        real(kind=realtype), dimension(3) :: mypt, otherpt, enorm
        real(kind=realtype) :: fact, dstar, curdist, mindist, edgelength
        integer(kind=intTYpe) :: otherID, otherIndex, closestOtherIndex, closestOtherString
        integer(kind=intType) :: id, index
        real(kind=realtype) :: timea, pt(3), v(3), costheta, cutoff, dist, maxh, ratio
 
        if (nstrings == 0) then
            return
        end if
 
        ! Now make we determine the nearest point on another substring
        ! for each point.
        nalloc = 50
        allocate (results(nalloc))
        master => strings(1)%p
 
        ! Loop over the fullStrings
        do i = 1, nstrings
            str => strings(i) ! Easier readability
 
            ! No need to do anything with the pocket string.
            if (str%isPocket) then
                cycle
            end if
 
            ! Allocate space for otherID as it is not done yet
            if (associated(str%otherID)) then
                deallocate (str%otherID)
            end if
 
            allocate (str%otherID(2, str%nNodes))
            str%otherID = -1
 
            ! Loop over my nodes and search for it in master tree
            nodeloop: do j = 1, str%nNodes
 
                ! Set the initial maximum number of neighbours
                ! This can be at most the total number of nodes
                nsearch = min(nalloc, master%nNodes)
 
                ! We have to be careful since single-sided chains have only
                ! 1 neighbour at each end.
 
                call getnodeinfo(str, j, checkleft, checkright, concave, &
                                 xj, xjm1, xjp1, normj)
                isendnode = .false.
                if (.not. (checkleft .eqv. checkright)) then
                    ! Since we don't need to check one side, this means we're
                    ! at the end of the chain. This is important since this
                    ! node *MUST* be attached to another node on another
                    ! chain at the end
                    isendnode = .true.
                end if
 
                outerloop: do
                    mindist = large
                    closestotherindex = -1
                    call kdtree2_n_nearest(master%tree, xj, nsearch, results)
 
                    ! Only check edges connected to nodes within the
                    ! distance the maximum element size of my self or the
                    ! closest node. We put in a fudge factor of 1.5.
 
                    innerloop: do k = 1, nsearch
 
                        ! Since we know the results are sorted, if the
                        ! distance(k) > than our current minDist, we can stop
                        ! since there is no possible way that any of the
                        ! remaining points can be closer given that the modified
                        ! D* is always larger than the original D
 
                        ! Extract current information to make things a little
                        ! easier to read
                        curdist = sqrt(results(k)%dis)
                        idx = results(k)%idx
                        pt = master%x(:, idx)
 
                        ! ---------------------------------------------
                        ! Exit Condition: We can stop the loop if the current
                        ! uncorrected distance is larger than our current
                        ! minimum. This guarantees the minimum corrected
                        ! distance is found.
                        ! ---------------------------------------------
 
                        if (curdist > mindist) then
                            exit outerloop
                        end if
 
                        ! ---------------------------------------------
                        ! Check 1: If the node we found isn't on our
                        ! substring. we don't need to do anything
                        ! ---------------------------------------------
 
                        if (master%cNodes(1, idx) == str%myID) then
                            cycle innerloop
                        end if
 
                        ! ---------------------------------------------
 
                        ! Check 1b: If the node we found has been removed due
                        ! to self zipping, we can just keep going
                        ! --------------------------------------------
                        if (master%cNodes(2, idx) == 0) then
                            cycle innerloop
                        end if
 
                        ! The first time we make it here, idx will be the
                        ! index of the closest node on another string that
                        ! isn't me.
                        if (closestotherindex == -1) then
                            closestotherstring = master%cNodes(1, idx)
                            closestotherindex = master%cNodes(2, idx)
                        end if
 
                        ! ---------------------------------------------
                        ! Check 2: Check if the node we found violates the
                        ! the "in front" test. For a concave corner TWO
                        ! triangle areas formed by the point and the two
                        ! edges must be positive. For a convex corner only
                        ! one of the triangle areas needs to be positive.
                        ! ---------------------------------------------
                        if (.not. nodeinfrontofedges(pt, concave, checkleft, checkright, &
                                                     xj, xjm1, xjp1, normj)) then
                            cycle innerloop
                        end if
 
                        ! ---------------------------------------------
                        ! Check 3: This is the *reverse* of check 2: Is the
                        ! node we're searching for visible from the potential
                        ! closest other node.
                        ! ---------------------------------------------
                        otherid = master%cNodes(1, idx)
                        otherindex = master%cNodes(2, idx)
 
                        call getnodeinfo(strings(otherid), otherindex, checkleft2, &
                                         checkright2, concave2, xk, xkm1, xkp1, normk)
 
                        if (.not. nodeinfrontofedges(xj, concave2, checkleft2, &
                                                     checkright2, xk, xkm1, xkp1, normk)) then
                            cycle innerloop
                        end if
 
                        ! ---------------------------------------------
                        ! Check 4a: Check if the potential node intersects
                        ! itself.
                        ! ---------------------------------------------
                        if (overlappededges(str, j, pt)) then
                            cycle
                        end if
 
                        ! ---------------------------------------------
                        ! Check 4b: OR if the other node would have to
                        ! intersect *ITSELF* to get back to me. This is used
                        ! to catch closest points crossing over thin strips.
                        ! ---------------------------------------------
 
                        if (overlappededges(strings(otherid), otherindex, xj)) then
                            cycle
                        end if
 
                        ! ---------------------------------------------
                        ! Check 4c: Make sure it doesn't inersect the closest
                        ! string if that happens to be different from the
                        ! cloest one.  string. This should only check very
                        ! rare cases the other checks miss.
                        ! ---------------------------------------------
 
                        if (otherid /= closestotherstring) then
                            if (overlappededges2( &
                                strings(closestotherstring), xj, normj, pt)) then
                                cycle
                            end if
                        end if
 
                        ! ---------------------------------------------
                        ! Check 4d: If this is an end node, we need to check
                        ! if the potential canditate is also a end node
                        ! ---------------------------------------------
                        if (isendnode) then
                            if (checkright2 .eqv. checkleft2) then
                                cycle
                            end if
                        end if
 
                        ! ---------------------------------------------
                        ! Check 5: Now that the point has passed the previous
                        ! checks, we can compute the agumented distance
                        ! function and see if it better than the exisitng min
                        ! distance.
                        ! ---------------------------------------------
 
                        ! Now calculate our new distance
                        v = pt - xj
                        v = v / mynorm2(v)
 
                        ! Recompute the distance function
                        costheta = abs(dot_product(normj, v))
 
                        ! Update distFunction
                        dstar = curdist / (max(1 - costheta, 1e-6))
 
                        if (dstar < mindist) then
                            ! Save the string ID and the index.
                            mindist = dstar
                            str%otherID(:, j) = master%cNodes(:, idx)
                        end if
                    end do innerloop
 
                    ! If we have already searched the max, we have to quit the loop
                    if (nsearch == master%Nnodes) then
                        exit outerloop
                    end if
 
                    ! We are not 100% sure that we found the minium
                    ! yet. Make nAlloc twice as big and start over.
                    nsearch = nsearch * 2
                    nsearch = min(nsearch, master%nNodes)
                    if (nsearch > nalloc) then
                        deallocate (results)
                        nalloc = nalloc * 2
                        allocate (results(nalloc))
                    end if
                end do outerloop
            end do nodeloop
 
            ! Do a sanity check to fix some extraordinary cases. If a node
            ! hasn't found a neighbouring string but each of the two nodes
            ! either side have, and they found the *same* string, just
            ! accept that.
 
            do j = 3, str%nNodes - 2
                if (str%otherID(1, j) == -1) then
                    ! Bad node:
                    oid(1) = str%otherID(1, j - 2)
                    oid(2) = str%otherID(1, j - 1)
                    oid(3) = str%otherID(1, j + 1)
                    oid(4) = str%otherID(1, j + 2)
 
                    if (oid(1) /= -1 .and. &
                        oid(1) == oid(2) .and. &
                        oid(1) == oid(3) .and. &
                        oid(1) == oid(4)) then
 
                        if (debugzipper) then
                            print *, '****************************************************************'
                            print *, 'Warning: Fixing a bad association on string ', i, 'at index', j
                            print *, '****************************************************************'
                        end if
 
                        ! We have a '-1' surrounded by the same gap string
 
                        ! Set the stringID
                        str%otherID(1, j) = oid(1)
 
                        ! Estimate what the other index should be. Since this
                        ! is in the middle of the string, the exact index
                        ! shouldn't matter.
                        str%otherID(2, j) = str%otherID(2, j - 1)
                    end if
                end if
            end do
 
        end do
    end subroutine stringmatch
 
    subroutine writeoversetstring(str, strings, n, fileID)
 
        use constants
        use utils, only: mynorm2
        use commonformats, only: sci12, int5
        implicit none
 
        integer(kind=intType), intent(in) :: fileID, n
        type(oversetstring), intent(in) :: str
        type(oversetstring), intent(in), dimension(n) :: strings
        integer(kind=intType) :: i, j, id, index
        real(kind=realtype), dimension(3) :: mypt, otherpt, vec
        real(kind=realtype) :: maxh, dist, ratio
 
        character(80) :: zoneName
 
        write (zonename, "(a,I5.5)") "Zone T=gap_", str%myID
        write (fileid, *) trim(zonename)
 
        write (fileid, *) "Nodes = ", str%nNodes, " Elements= ", str%nElems, " ZONETYPE=FELINESEG"
        write (fileid, *) "DATAPACKING=BLOCK"
 
        ! Nodes
        do j = 1, 3
            do i = 1, str%nNodes
                write (fileid, sci12) str%x(j, i)
            end do
        end do
 
        ! Node normal
        do j = 1, 3
            do i = 1, str%nNodes
                write (fileid, sci12) str%norm(j, i)
            end do
        end do
 
        ! Vector between closest points
        do j = 1, 3
            do i = 1, str%nNodes
                mypt = str%x(:, i)
                id = str%otherID(1, i)
                if (id /= -1) then
                    index = str%otherID(2, i)
                    otherpt = strings(id)%x(:, index)
                    vec = otherpt - mypt
                else
                    vec = zero
                end if
 
                write (fileid, sci12) vec(j)
            end do
        end do
 
        ! global node ID
        do i = 1, str%nNodes
            write (fileid, sci12) real(str%ind(i))
        end do
 
        ! gapID
        do i = 1, str%nNodes
            write (fileid, sci12) real(str%myID)
        end do
 
        ! gap Index
        do i = 1, str%nNodes
            write (fileid, sci12) real(i)
        end do
 
        if (associated(str%otherID)) then
            ! otherID
            do i = 1, str%nNodes
                write (fileid, sci12) real(str%otherID(1, i))
            end do
 
            ! other Index
            do i = 1, str%nNodes
                write (fileid, sci12) real(str%otherID(2, i))
            end do
        else
            do i = 1, 2 * str%nNodes
                write (fileid, sci12) zero
            end do
        end if
 
        do i = 1, str%nNodes
            mypt = str%x(:, i)
            id = str%otherID(1, i)
            if (id /= -1) then
                index = str%otherID(2, i)
                otherpt = strings(id)%x(:, index)
                dist = mynorm2(mypt - otherpt)
                maxh = max(str%h(i), strings(id)%h(index))
                ratio = dist / maxh
            else
                ratio = zero
            end if
 
            write (fileid, sci12) ratio
        end do
 
        do i = 1, str%nElems
            write (fileid, int5) str%conn(1, i), str%conn(2, i)
        end do
 
    end subroutine writeoversetstring
 
    subroutine writeoversetmaster(str, fileID)
 
        use constants
        use utils, only: mynorm2
        use commonformats, only: sci12, int5
        implicit none
 
        type(oversetstring), intent(in) :: str
        integer(kind=intType), intent(in) :: fileID
        integer(kind=intType) :: i, j, id, index
        real(kind=realtype), dimension(3) :: mypt, otherpt, vec
        real(kind=realtype) :: maxh, dist, ratio
 
        character(80) :: zoneName
 
        write (zonename, "(a,I5.5)") "Zone T=gap_", str%myID
        write (fileid, *) trim(zonename)
 
        write (fileid, *) "Nodes = ", str%nNodes, " Elements= ", str%nElems, " ZONETYPE=FELINESEG"
        write (fileid, *) "DATAPACKING=BLOCK"
 
        ! Nodes
        do j = 1, 3
            do i = 1, str%nNodes
                write (fileid, sci12) str%x(j, i)
            end do
        end do
 
        do i = 1, str%nElems
            write (fileid, int5) str%conn(1, i), str%conn(2, i)
        end do
 
    end subroutine writeoversetmaster
 
    subroutine writeoversettriangles(string, fileName, startTri, endTri)
 
        use constants
        use commonformats, only: sci12
        implicit none
 
        type(oversetstring), intent(inout) :: string
        integer(kind=intType), intent(in) :: startTri, endTri
        character(*) :: fileName
        integer(kind=intType) :: i, j
        character(80) :: zoneName
 
        open (unit=101, file=trim(filename), form='formatted')
        write (101, *) 'TITLE = "Triangles"'
        write (101, *) 'Variables = "X", "Y", "Z"'
 
        write (zonename, "(a,I5.5)") "Zone T=triangles_", string%myID
        write (101, *) trim(zonename)
 
        write (101, *) "Nodes = ", string%nNodes, " Elements= ", (endtri - starttri + 1), " ZONETYPE=FETRIANGLE"
        write (101, *) "DATAPACKING=POINT"
 
        ! Write all the coordinates
        do i = 1, string%nNodes
            do j = 1, 3
                write (101, sci12, advance='no') string%x(j, i)
            end do
            write (101, "(1x)")
        end do
 
        do i = starttri, endtri
            write (101, "(*(I7))") string%tris(1, i), string%tris(2, i), string%tris(3, i)
        end do
        close (101)
    end subroutine writeoversettriangles
 
    subroutine writezipperdebug(str)
 
        ! Save the state of an unsplit string such that it can be debugged
        ! later without running overset interpolation.
 
        use constants
        implicit none
 
        type(oversetstring) :: str
        integer(kind=intType) :: i, j
 
        open (unit=101, file="debug.zipper", form='formatted')
        write (101, *) str%nNodes
        write (101, *) str%nElems
        do i = 1, str%nNodes
            do j = 1, 10
                write (101, *) str%nodeData(j, i)
            end do
        end do
 
        do i = 1, str%nElems
            do j = 1, 2
                write (101, *) str%conn(j, i)
            end do
        end do
 
        do i = 1, str%nNodes
            do j = 1, 3
                write (101, *) str%intNodeData(j, i)
            end do
        end do
        close (101)
    end subroutine writezipperdebug
 
    subroutine loadzipperdebug(fileName, str)
 
        ! Save the state of an unsplit string such that it can be debugged
        ! later without running overset interpolation.
 
        use constants
        implicit none
 
        character(*), intent(in) :: fileName
        type(oversetstring) :: str
        integer(kind=intType) :: i, j
 
        open (unit=101, file=filename, form='formatted')
        read (101, *) str%nNodes
        read (101, *) str%nElems
        call nullifystring(str)
 
        allocate (str%nodeData(10, str%nNodes))
        allocate (str%conn(2, str%nElems))
        allocate (str%intNodeData(3, str%nNodes))
 
        do i = 1, str%nNodes
            do j = 1, 10
                read (101, *) str%nodeData(j, i)
            end do
        end do
 
        do i = 1, str%nElems
            do j = 1, 2
                read (101, *) str%conn(j, i)
            end do
        end do
 
        do i = 1, str%nNodes
            do j = 1, 3
                read (101, *) str%intNodeData(j, i)
            end do
        end do
        close (101)
 
        call setstringpointers(str)
 
    end subroutine loadzipperdebug
 
    subroutine pointintriangle(x1, x2, x3, pt, inTri)
 
        use constants
        use utils, only: cross_prod
        implicit none
        real(kind=realtype), dimension(3), intent(in) :: x1, x2, x3, pt
        logical, intent(out) :: inTri
 
        if (sameside(pt, x1, x2, x3) .and. sameside(pt, x2, x1, x3) .and. sameside(pt, x3, x1, x2)) then
            intri = .true.
        else
            intri = .false.
        end if
 
    contains
        function sameside(p1, p2, a, b)
 
            implicit none
            logical :: sameside
            real(kind=realtype), dimension(3) :: p1, p2, a, b, cp1, cp2
 
            sameside = .false.
            call cross_prod(b - a, p1 - a, cp1)
            call cross_prod(b - a, p2 - a, cp2)
            if (dot_product(cp1, cp2) >= zero) then
                sameside = .true.
            end if
        end function sameside
    end subroutine pointintriangle
 
    function positivetriarea(p1, p2, p3, norm)
 
        use constants
        use utils, only: cross_prod
        implicit none
        real(kind=realtype), intent(in), dimension(3) :: p1, p2, p3, norm
        real(kind=realtype), dimension(3) :: n
        logical :: positivetriarea
 
        call cross_prod(p2 - p1, p3 - p1, n)
        if (dot_product(n, norm) > zero) then
            positivetriarea = .true.
        else
            positivetriarea = .false.
        end if
    end function positivetriarea
 
    subroutine getnodeinfo(str, j, checkLeft, checkRight, concave, xj, xjm1, xjp1, normj)
 
        use constants
        use utils, only: cross_prod
        implicit none
 
        type(oversetstring) :: str
        integer(kind=intType) :: j
        logical :: checkLeft, checkRight, concave
        real(kind=realtype), dimension(3) :: xj, xjm1, xjp1, normj
        real(kind=realtype), dimension(3) :: v
        checkleft = .true.
        checkright = .true.
        xj = str%x(:, j)
        normj = str%norm(:, j)
        concave = .false.
        if (str%isPeriodic) then
            if (j > 1 .and. j < str%nNodes) then
                xjm1 = str%x(:, j - 1)
                xjp1 = str%x(:, j + 1)
            else if (j == 1) then
                xjm1 = str%x(:, str%nNodes)
                xjp1 = str%x(:, j + 1)
            else if (j == str%nNodes) then
                xjm1 = str%x(:, j - 1)
                xjp1 = str%x(:, 1)
            end if
        else
            ! Not periodic. Assume the ends are concave. This will
            ! forces checking if both the left and right are ok,
            ! which since the leftOK and rightOK's default to
            ! .True., it just checks the one triangle which is what
            ! we want.
            if (j == 1) then
                checkleft = .false.
                concave = .true.
            end if
 
            if (j == str%nNodes) then
                checkright = .false.
                concave = .true.
            end if
 
            if (checkleft) &
                xjm1 = str%x(:, j - 1)
            if (checkright) &
                xjp1 = str%x(:, j + 1)
        end if
 
        if (checkleft .and. checkright) then
 
            ! Determine if the point is convex or concave provided
            !  we have both neighbours.
            call cross_prod(xjm1 - xj, xjp1 - xj, v)
 
            if (dot_product(v, normj) > zero) then
                concave = .true.
            end if
        end if
 
    end subroutine getnodeinfo
 
    function nodeinfrontofedges(pt, concave, checkLeft, checkRight, xj, xjm1, xjp1, normj)
 
        use constants
        implicit none
 
        real(kind=realtype), dimension(3), intent(in) :: pt, xj, xjm1, xjp1, normj
        logical, intent(in) :: concave, checkleft, checkright
        logical :: nodeinfrontofedges
        logical :: leftok, rightok
 
        nodeinfrontofedges = .true.
        leftok = .true.
        rightok = .true.
        if (checkleft .and. .not. positivetriarea(xj, xjm1, pt, normj)) then
            leftok = .false.
        end if
 
        if (checkright .and. .not. positivetriarea(xjp1, xj, pt, normj)) then
            rightok = .false.
        end if
 
        if (concave) then
            if (.not. (leftok .and. rightok)) then
                nodeinfrontofedges = .false.
            end if
        else
            if (.not. (leftok .or. rightok)) then
                nodeinfrontofedges = .false.
            end if
        end if
    end function nodeinfrontofedges
 
    function overlappededges(str, j, pt)
 
        use constants
        use utils, only: mynorm2, cross_prod
 
        implicit none
 
        ! Input/output
        real(kind=realtype), dimension(3), intent(in) :: pt
        type(oversetstring), intent(in) :: str
        integer(kind=intType), intent(in) :: j
        logical :: overlappededges
 
        ! Working
        integer(kind=intType) :: i
        real(kind=realtype), dimension(3) :: v, p1, p2, u, norma, normb, x0, norm
        real(kind=realtype) :: unrm, x1, x2, x3, x4, y1, y2, y3, y4, idet, px, py
        real(kind=realtype) :: u1, u2, v1, v2, w1, w2
        real(kind=realtype) :: s1, s2, tmp, line(2), vec(2), tol
        overlappededges = .false.
        tol = 1e-6
        ! We will conver this completely into a 2D problem by projecting
        ! everything onto the plane defined by norm. x0 is at the origin of
        ! the 2D system and the xaxis point from x0 to pt
 
        x0 = str%x(:, j)
        norm = str%norm(:, j)
 
        u = pt - x0
        unrm = mynorm2(u)
        u = u / unrm
 
        call cross_prod(norm, u, v)
        v = v / mynorm2(v)
 
        ! Now u,v,norm is an orthogonal coordinate system
        x1 = zero
        y1 = zero
        x2 = unrm
        y2 = zero
        overlappededges = .false.
 
        ! Loop over the number of edges on my string
        elemloop: do i = 1, str%nElems
            ! Don't check the ones right next to me, since they will
            ! "overlap" exactly at x0
 
            if (str%conn(1, i) == j .or. str%conn(2, i) == j) then
                cycle
            end if
 
            ! Project the two points into the plane
            p1 = str%x(:, str%conn(1, i))
            p2 = str%x(:, str%conn(2, i))
 
            norma = str%norm(:, str%conn(1, i))
            normb = str%norm(:, str%conn(2, i))
 
            ! Make sure the edges are on the same plane, otherwise this is
            ! meaningless
            if (dot_product(norma, norm) < half .or. dot_product(normb, norm) < half) then
                cycle
            end if
            ! Project the two points onto the plane
            p1 = p1 - norm * dot_product(p1 - x0, norm)
            p2 = p2 - norm * dot_product(p2 - x0, norm)
 
            ! Now get the 2D coordinates
            x3 = dot_product(p1 - x0, u)
            y3 = dot_product(p1 - x0, v)
            x4 = dot_product(p2 - x0, u)
            y4 = dot_product(p2 - x0, v)
 
            u1 = x2 - x1
            y2 = y2 - y1
 
            v1 = x4 - x3
            v2 = y4 - y3
 
            w1 = x1 - x3
            w2 = y1 - y3
 
            s1 = (v2 * w1 - v1 * w2) / (v1 * u2 - v2 * u1)
            s2 = (u1 * w2 - u2 * w1) / (u1 * v2 - u2 * v1)
 
            if (s1 > tol .and. s1 < one - tol .and. s2 > tol .and. s2 < one - tol) then
                overlappededges = .true.
                exit elemloop
            end if
        end do elemloop
 
    end function overlappededges
 
    function overlappededges2(str, pt1, norm, pt2)
 
        use constants
        use utils, only: mynorm2, cross_prod
        implicit none
 
        ! Input/output
        real(kind=realtype), dimension(3), intent(in) :: pt1, pt2, norm
        type(oversetstring), intent(in) :: str
        logical :: overlappededges2
 
        ! Working
        integer(kind=intType) :: i
        real(kind=realtype), dimension(3) :: v, p1, p2, u, norma, normb, x0
        real(kind=realtype) :: unrm, x1, x2, x3, x4, y1, y2, y3, y4, idet, px, py
        real(kind=realtype) :: u1, u2, v1, v2, w1, w2
        real(kind=realtype) :: s1, s2, tmp, line(2), vec(2), tol
 
        tol = 1e-6
        ! We will conver this completely into a 2D problem by projecting
        ! everything onto the plane defined by norm. x0 is at the origin of
        ! the 2D system and the xaxis point from x0 to pt
 
        x0 = pt1
 
        u = pt2 - x0
        unrm = mynorm2(u)
        u = u / unrm
 
        call cross_prod(norm, u, v)
        v = v / mynorm2(v)
 
        ! Now u,v,norm is an orthogonal coordinate system
        x1 = zero
        y1 = zero
        x2 = unrm
        y2 = zero
        overlappededges2 = .false.
 
        ! Loop over the number of edges on my string
        elemloop: do i = 1, str%nElems
 
            ! Project the two points into the plane
            p1 = str%x(:, str%conn(1, i))
            p2 = str%x(:, str%conn(2, i))
 
            norma = str%norm(:, str%conn(1, i))
            normb = str%norm(:, str%conn(2, i))
 
            ! Make sure the edges are on the same plane, otherwise this is
            ! meaningless
            if (dot_product(norma, norm) < half .or. dot_product(normb, norm) < half) then
                cycle
            end if
            ! Project the two points onto the plane
            p1 = p1 - norm * dot_product(p1 - x0, norm)
            p2 = p2 - norm * dot_product(p2 - x0, norm)
 
            ! Now get the 2D coordinates
            x3 = dot_product(p1 - x0, u)
            y3 = dot_product(p1 - x0, v)
            x4 = dot_product(p2 - x0, u)
            y4 = dot_product(p2 - x0, v)
 
            u1 = x2 - x1
            y2 = y2 - y1
 
            v1 = x4 - x3
            v2 = y4 - y3
 
            w1 = x1 - x3
            w2 = y1 - y3
 
            s1 = (v2 * w1 - v1 * w2) / (v1 * u2 - v2 * u1)
            s2 = (u1 * w2 - u2 * w1) / (u1 * v2 - u2 * v1)
 
            if (s1 > tol .and. s1 < one - tol .and. s2 > tol .and. s2 < one - tol) then
                overlappededges2 = .true.
                exit elemloop
            end if
        end do elemloop
 
    end function overlappededges2
 
end module stringops