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
    integer :: iEl,j,e2,bcid
    integer :: count,n
    integer,allocatable :: elems(:),sds(:)

    ! 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(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(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)
      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(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(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)
      endif
      bc => bc%next
    enddo

  endsubroutine MapBoundaryConditions_DGModel2D_t