RefineConnectivity Subroutine

  subroutine RefineConnectivity(nElem,baseSideInfo,baseCorner,nodeOffset,nUniqueSidesBase, &
                                refSideInfo,refCorner,nUniqueSidesRef)
    !! Build the connectivity of a uniformly (2:1 in each direction) refined 2-D mesh from its
    !! base connectivity. Every base element p is split into four children (global id
    !! 4*(p-1)+c). The scheme is fully deterministic integer bookkeeping - no coordinate
    !! hashing, no flip recomputation, no MPI:
    !!
    !!   - Sibling faces interior to a parent are same-orientation (flip 0).
    !!   - A child face on a parent boundary inherits the parent face's neighbor element (the
    !!     appropriate child of the parent's neighbor), neighbor side, and flip. The sub-position
    !!     pairing across the face uses the parent flip: sub-position t maps to t on the neighbor
    !!     when flip = 0 and to 3-t when flip = 1.
    !!
    !! Because refinement is orientation preserving, this reproduces exactly the flips that a
    !! corner-node matching pass (RecalculateFlip) would compute, but without needing globally
    !! consistent node ids, which structured meshes do not guarantee.
    !!
    !! Refined corner node ids (returned in refCorner, SW/SE/NE/NW per child) reuse the base
    !! corner ids for retained corners, place one shared midpoint per base side
    !! (nodeOffset + base global side id) and one center per base element
    !! (nodeOffset + nUniqueSidesBase + p). These feed globalNodeIDs for mesh I/O; the solver
    !! reads flips from refSideInfo directly.
    implicit none
    integer,intent(in) :: nElem
    integer,intent(in) :: baseSideInfo(1:5,1:4,1:nElem)
    integer,intent(in) :: baseCorner(1:4,1:nElem)
    integer,intent(in) :: nodeOffset
    integer,intent(in) :: nUniqueSidesBase
    integer,intent(out) :: refSideInfo(1:5,1:4,1:4*nElem)
    integer,intent(out) :: refCorner(1:4,1:4*nElem)
    integer,intent(out) :: nUniqueSidesRef
    ! Local
    ! Interior (sibling) face pairs within a parent: (childA,sideA) <-> (childB,sideB), flip 0.
    integer,parameter :: nInner = 4
    integer,parameter :: innerA(1:2,1:nInner) = reshape([1,2,1,3,2,3,4,2],[2,nInner])
    integer,parameter :: innerB(1:2,1:nInner) = reshape([2,4,4,1,3,1,3,4],[2,nInner])
    integer :: p,s,t,tq,k,e,nbr,nbrSide
    integer :: q,s2,f,childP,childQ,eBase
    integer :: m(1:4),z,cc(1:4)

    refSideInfo = 0
    refCorner = 0

    do p = 1,nElem
      eBase = 4*(p-1)

      ! ---- Corner node ids for the four children ----
      cc(1:4) = baseCorner(1:4,p)
      do s = 1,4
        m(s) = nodeOffset+abs(baseSideInfo(2,s,p)) ! shared midpoint on base side s
      enddo
      z = nodeOffset+nUniqueSidesBase+p ! element center
      refCorner(1:4,eBase+1) = [cc(1),m(1),z,m(4)] ! child 1 (SW)
      refCorner(1:4,eBase+2) = [m(1),cc(2),m(2),z] ! child 2 (SE)
      refCorner(1:4,eBase+3) = [z,m(2),cc(3),m(3)] ! child 3 (NE)
      refCorner(1:4,eBase+4) = [m(4),z,m(3),cc(4)] ! child 4 (NW)

      ! ---- Interior sibling faces (flip 0) ----
      do k = 1,nInner
        childP = innerA(1,k); s = innerA(2,k)
        refSideInfo(3,s,eBase+childP) = eBase+innerB(1,k) ! neighbor element
        refSideInfo(4,s,eBase+childP) = 10*innerB(2,k) ! 10*neighbor side + flip(0)
        refSideInfo(5,s,eBase+childP) = 0 ! interior: no BC
        childP = innerB(1,k); s = innerB(2,k)
        refSideInfo(3,s,eBase+childP) = eBase+innerA(1,k)
        refSideInfo(4,s,eBase+childP) = 10*innerA(2,k)
        refSideInfo(5,s,eBase+childP) = 0
      enddo

      ! ---- Exterior faces on each parent side ----
      do s = 1,4
        q = baseSideInfo(3,s,p) ! neighbor base element (0 = physical boundary)
        s2 = baseSideInfo(4,s,p)/10 ! neighbor base side
        f = baseSideInfo(4,s,p)-10*s2 ! base flip (0 or 1)
        do t = 1,2
          childP = childOfSide(t,s) ! child on this side; its local side is s
          if(q == 0) then
            ! Physical boundary: inherit the parent side's BC, no neighbor.
            refSideInfo(3,s,eBase+childP) = 0
            refSideInfo(4,s,eBase+childP) = 0
            refSideInfo(5,s,eBase+childP) = baseSideInfo(5,s,p)
          else
            if(f == 0) then
              tq = t
            else
              tq = 3-t
            endif
            childQ = childOfSide(tq,s2)
            refSideInfo(3,s,eBase+childP) = 4*(q-1)+childQ
            refSideInfo(4,s,eBase+childP) = 10*s2+f
            refSideInfo(5,s,eBase+childP) = 0
          endif
        enddo
      enddo
    enddo

    ! ---- Assign a global side id to each unique refined side (count each shared side once) ----
    nUniqueSidesRef = 0
    do e = 1,4*nElem
      do s = 1,4
        if(refSideInfo(2,s,e) /= 0) cycle ! already numbered from its partner
        nbr = refSideInfo(3,s,e)
        nbrSide = refSideInfo(4,s,e)/10
        nUniqueSidesRef = nUniqueSidesRef+1
        refSideInfo(2,s,e) = nUniqueSidesRef
        if(nbr /= 0) refSideInfo(2,nbrSide,nbr) = nUniqueSidesRef
      enddo
    enddo

    ! Side type field is unused in SELF's 2-D path; leave refSideInfo(1,:,:) = 0.

  endsubroutine RefineConnectivity