MapBoundaryConditions_DGModel2D_t Subroutine

  subroutine MapBoundaryConditions_DGModel2D_t(this)
    !! Scan the mesh sideInfo and populate the elements/sides
    !! arrays for each registered boundary condition.
    implicit none
    class(DGModel2D_t),intent(inout) :: this
    ! Local
    type(BoundaryCondition),pointer :: bc
    type(BoundaryCondition),pointer :: bcnode
    integer :: iEl,j,e2,bcid
    integer :: count,n
    integer :: nUnmapped,idUnmapped,iError
    integer :: idRank,srcRank
    logical :: skipMortars
    integer,allocatable :: elems(:),sds(:)

    ! Mortar sides carry sideInfo(3) = 0 but are interior faces; sideInfo(1) holds the mortar
    ! index (see SELF_Mesh_2D_t / SELF_Mesh_3D_t). They must be excluded from every pass below,
    ! forward and reverse alike: a model may legitimately register bcid 0, and a mortar side
    ! carries sideInfo(5) = 0, so a forward pass matching on sideInfo(5) alone would put
    ! interior faces into that condition's list and SetBoundaryCondition would then overwrite
    ! what the mortar exchange had just written. The HOPr reader copies sideInfo(1) verbatim
    ! from the file, where it is the HOPr side type and may be nonzero on an ordinary face, so
    ! sideInfo(1) is only a mortar marker on a mesh that actually carries mortars.
    skipMortars = this%mesh%nMortars > 0

    ! Map hyperbolic BCs
    bc => this%hyperbolicBCs%head
    do while(associated(bc))
      ! Pass 1: count boundary faces for this bcid
      count = 0
      do iEl = 1,this%mesh%nElem
        do j = 1,4
          e2 = this%mesh%sideInfo(3,j,iEl)
          bcid = this%mesh%sideInfo(5,j,iEl)
          if(skipMortars .and. this%mesh%sideInfo(1,j,iEl) /= 0) cycle
          if(e2 == 0 .and. bcid == bc%bcid) count = count+1
        enddo
      enddo

      if(count > 0) then
        ! Pass 2: fill element/side arrays
        allocate(elems(count),sds(count))
        n = 0
        do iEl = 1,this%mesh%nElem
          do j = 1,4
            e2 = this%mesh%sideInfo(3,j,iEl)
            bcid = this%mesh%sideInfo(5,j,iEl)
            if(skipMortars .and. this%mesh%sideInfo(1,j,iEl) /= 0) cycle
            if(e2 == 0 .and. bcid == bc%bcid) then
              n = n+1
              elems(n) = iEl
              sds(n) = j
            endif
          enddo
        enddo
        call this%hyperbolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds)
        deallocate(elems,sds)
      else
        ! Drop any mapping left by a previous call. SetBoundaryCondition dispatches every
        ! registered condition and each one loops over its own nBoundaries, so a stale
        ! element/side list would keep this condition writing faces it no longer owns once
        ! the mesh is re-tagged - and the unmapped tally below would then describe something
        ! other than what runs. Clearing nBoundaries also parks the GPU wrappers, which guard
        ! on the same field, so the device arrays are never read while stale.
        bc%nBoundaries = 0
        if(allocated(bc%elements)) deallocate(bc%elements)
        if(allocated(bc%sides)) deallocate(bc%sides)
      endif
      bc => bc%next
    enddo

    ! Map parabolic BCs
    bc => this%parabolicBCs%head
    do while(associated(bc))
      count = 0
      do iEl = 1,this%mesh%nElem
        do j = 1,4
          e2 = this%mesh%sideInfo(3,j,iEl)
          bcid = this%mesh%sideInfo(5,j,iEl)
          if(skipMortars .and. this%mesh%sideInfo(1,j,iEl) /= 0) cycle
          if(e2 == 0 .and. bcid == bc%bcid) count = count+1
        enddo
      enddo

      if(count > 0) then
        allocate(elems(count),sds(count))
        n = 0
        do iEl = 1,this%mesh%nElem
          do j = 1,4
            e2 = this%mesh%sideInfo(3,j,iEl)
            bcid = this%mesh%sideInfo(5,j,iEl)
            if(skipMortars .and. this%mesh%sideInfo(1,j,iEl) /= 0) cycle
            if(e2 == 0 .and. bcid == bc%bcid) then
              n = n+1
              elems(n) = iEl
              sds(n) = j
            endif
          enddo
        enddo
        call this%parabolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds)
        deallocate(elems,sds)
      else
        ! Drop any mapping left by a previous call. SetBoundaryCondition dispatches every
        ! registered condition and each one loops over its own nBoundaries, so a stale
        ! element/side list would keep this condition writing faces it no longer owns once
        ! the mesh is re-tagged - and the unmapped tally below would then describe something
        ! other than what runs. Clearing nBoundaries also parks the GPU wrappers, which guard
        ! on the same field, so the device arrays are never read while stale.
        bc%nBoundaries = 0
        if(allocated(bc%elements)) deallocate(bc%elements)
        if(allocated(bc%sides)) deallocate(bc%sides)
      endif
      bc => bc%next
    enddo

    ! Reverse check. Both loops above iterate over registrations, so a boundary edge whose
    ! bcid matches no registration is never enumerated and its exterior state is never
    ! written. Count those edges here: the sideInfo scan is already what this routine costs,
    ! and it runs at Init and after every Regrid, never inside the time loop.
    !
    nUnmapped = 0
    idUnmapped = -1
    do iEl = 1,this%mesh%nElem
      do j = 1,4
        e2 = this%mesh%sideInfo(3,j,iEl)
        if(e2 /= 0) cycle ! interior or rank-shared: sideInfo(3) holds a global element id
        if(skipMortars .and. this%mesh%sideInfo(1,j,iEl) /= 0) cycle
        bcid = this%mesh%sideInfo(5,j,iEl)
        ! Only the HYPERBOLIC list decides whether the face is handled. SetBoundaryCondition
        ! dispatches that list alone and it is what writes solution%extBoundary, the trace the
        ! Riemann solver consumes; the parabolic list writes solutionGradient%extBoundary
        ! through SetGradientBoundaryCondition. A bcid registered only parabolically therefore
        ! leaves the solution trace unwritten - exactly the failure this scan exists to catch.
        bcnode => this%hyperbolicBCs%GetBCForID(bcid)
        if(associated(bcnode)) cycle
        nUnmapped = nUnmapped+1
        ! Keep the FIRST offender, not the largest: a bcid is any integer, so a max()
        ! against a sentinel would never report one that sits below the sentinel.
        if(nUnmapped == 1) idUnmapped = bcid
      enddo
    enddo

    ! Each rank owns a slice of the mesh, so a bcid that appears nowhere here may still be
    ! present on another rank. Every rank reaches this routine on both the Init and the Regrid
    ! path, so the collective is safe.
    if(this%mesh%decomp%mpiEnabled) then
      call mpi_allreduce(nUnmapped,this%nUnmappedBoundaries,1,MPI_INTEGER, &
                         MPI_SUM,this%mesh%decomp%mpiComm,iError)
      ! A bcid is any integer, so no value can serve as an "absent" sentinel in a reduction
      ! over the id itself. Agree on the lowest-numbered rank that actually has an offender
      ! and take its id from there; a rank with none bids nRanks and so never wins.
      if(nUnmapped > 0) then
        idRank = this%mesh%decomp%rankId
      else
        idRank = this%mesh%decomp%nRanks
      endif
      call mpi_allreduce(idRank,srcRank,1,MPI_INTEGER,MPI_MIN, &
                         this%mesh%decomp%mpiComm,iError)
      if(srcRank < this%mesh%decomp%nRanks) then
        this%unmappedBoundaryID = idUnmapped
        call mpi_bcast(this%unmappedBoundaryID,1,MPI_INTEGER,srcRank, &
                       this%mesh%decomp%mpiComm,iError)
      else
        this%unmappedBoundaryID = -1
      endif
    else
      this%nUnmappedBoundaries = nUnmapped
      this%unmappedBoundaryID = idUnmapped
    endif
    ! Regrid remaps onto a new mesh, so a mesh that is still mis-tagged warns again.
    this%unmappedBoundariesReported = .false.

  endsubroutine MapBoundaryConditions_DGModel2D_t