MapBoundaryConditions_DGModel1D_t Subroutine

  subroutine MapBoundaryConditions_DGModel1D_t(this)
    !! Scan the mesh boundary condition IDs and populate the elements/sides
    !! arrays for each registered boundary condition.
    implicit none
    class(DGModel1D_t),intent(inout) :: this
    ! Local
    type(BoundaryCondition),pointer :: bc
    type(BoundaryCondition),pointer :: bcnode
    integer :: nelem,count,n
    integer :: endpoint,bcid
    integer :: elems(2),sds(2)

    nelem = this%mesh%nElem

    ! Map hyperbolic BCs
    bc => this%hyperbolicBCs%head
    do while(associated(bc))
      count = 0
      if(this%mesh%bcid(1) == bc%bcid) count = count+1
      if(this%mesh%bcid(2) == bc%bcid) count = count+1

      if(count > 0) then
        n = 0
        if(this%mesh%bcid(1) == bc%bcid) then
          n = n+1
          elems(n) = 1
          sds(n) = 1
        endif
        if(this%mesh%bcid(2) == bc%bcid) then
          n = n+1
          elems(n) = nelem
          sds(n) = 2
        endif
        call this%hyperbolicBCs%PopulateBoundaries(bc%bcid,count, &
                                                   elems(1:count),sds(1:count))
      else
        ! Drop any mapping left by a previous call. SetBoundaryCondition dispatches every
        ! registered condition over its own nBoundaries, so a stale element/side list would
        ! keep this condition writing endpoints it no longer owns once the mesh is re-tagged,
        ! and the tally below would then describe something other than what runs.
        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
      if(this%mesh%bcid(1) == bc%bcid) count = count+1
      if(this%mesh%bcid(2) == bc%bcid) count = count+1

      if(count > 0) then
        n = 0
        if(this%mesh%bcid(1) == bc%bcid) then
          n = n+1
          elems(n) = 1
          sds(n) = 1
        endif
        if(this%mesh%bcid(2) == bc%bcid) then
          n = n+1
          elems(n) = nelem
          sds(n) = 2
        endif
        call this%parabolicBCs%PopulateBoundaries(bc%bcid,count, &
                                                  elems(1:count),sds(1:count))
      else
        ! Drop any mapping left by a previous call. SetBoundaryCondition dispatches every
        ! registered condition over its own nBoundaries, so a stale element/side list would
        ! keep this condition writing endpoints it no longer owns once the mesh is re-tagged,
        ! and the tally below would then describe something other than what runs.
        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 an endpoint whose bcid
    ! matches no registration is never enumerated. A bcid of 0 is the deliberate periodic
    ! default (Mesh1D initialises bcid to 0), so it is not a fault.
    !
    !
    ! Mesh1D is replicated on every rank (nGlobalElem = nElem, no element decomposition), so
    ! the count is already global and no reduction is needed here.
    this%nUnmappedBoundaries = 0
    this%unmappedBoundaryID = -1
    do endpoint = 1,2
      bcid = this%mesh%bcid(endpoint)
      if(bcid == 0) cycle ! periodic by default
      ! Only the HYPERBOLIC list decides whether the endpoint is handled. SetBoundaryCondition
      ! dispatches that list alone and it is what writes solution%extBoundary; the parabolic
      ! list writes solutionGradient%extBoundary through SetGradientBoundaryCondition. A bcid
      ! registered only parabolically therefore leaves the solution trace on the periodic
      ! default - exactly the failure this scan exists to catch.
      bcnode => this%hyperbolicBCs%GetBCForID(bcid)
      if(associated(bcnode)) cycle
      this%nUnmappedBoundaries = this%nUnmappedBoundaries+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(this%nUnmappedBoundaries == 1) this%unmappedBoundaryID = bcid
    enddo
    this%unmappedBoundariesReported = .false.

  endsubroutine MapBoundaryConditions_DGModel1D_t