RefineConnectivity Subroutine

  subroutine RefineConnectivity(nElem,baseSideInfo,refSideInfo,nUniqueSidesRef)
    !! Build the connectivity of a uniformly (2:1 in each direction) refined 3-D mesh from
    !! its base connectivity. Every base element p is split into eight children (global id
    !! 8*(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 face, and flip. The
    !!     quadrant pairing across the face applies the parent flip's quadrant permutation
    !!     (faceQuadPerm).
    !!
    !! 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
    !! globalNodeIDs are not produced: the solver reads flips from refSideInfo directly and
    !! no 3-D consumer requires them (the adaptive-mesh emitter sets them to zero as well).
    implicit none
    integer,intent(in) :: nElem
    integer,intent(in) :: baseSideInfo(1:5,1:6,1:nElem)
    integer,intent(out) :: refSideInfo(1:5,1:6,1:8*nElem)
    integer,intent(out) :: nUniqueSidesRef
    ! Local
    ! Interior (sibling) face pairs within a parent: (childA,faceA) <-> (childB,faceB),
    ! flip 0. Four pairs per direction : East/West, North/South, Top/Bottom.
    integer,parameter :: nInner = 12
    integer,parameter :: innerA(1:2,1:nInner) = reshape([ &
                                                        1,3,4,3,5,3,8,3, & ! East faces
                                                        1,4,2,4,5,4,6,4, & ! North faces
                                                        1,6,2,6,3,6,4,6],[2,nInner]) ! Top faces
    integer,parameter :: innerB(1:2,1:nInner) = reshape([ &
                                                        2,5,3,5,6,5,7,5, & ! West faces
                                                        4,2,3,2,8,2,7,2, & ! South faces
                                                        5,1,6,1,7,1,8,1],[2,nInner]) ! Bottom faces
    integer :: p,s,t,tq,k,e,nbr,nbrSide
    integer :: q,s2,f,childP,childQ,eBase

    refSideInfo = 0

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

      ! ---- 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 face + 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 face ----
      do s = 1,6
        q = baseSideInfo(3,s,p) ! neighbor base element (0 = physical boundary)
        s2 = baseSideInfo(4,s,p)/10 ! neighbor base face
        f = baseSideInfo(4,s,p)-10*s2 ! base flip (0..7)
        do t = 1,4
          childP = childOfFace(t,s) ! child on this face; its local face is s
          if(q == 0) then
            ! Physical boundary: inherit the parent face'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
            tq = faceQuadPerm(t,f)
            childQ = childOfFace(tq,s2)
            refSideInfo(3,s,eBase+childP) = 8*(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 face (count each shared face once) ----
    nUniqueSidesRef = 0
    do e = 1,8*nElem
      do s = 1,6
        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 only used to mark mortar faces; uniform refinement of a
    ! conforming mesh is conforming, so leave refSideInfo(1,:,:) = 0.

  endsubroutine RefineConnectivity