load_BC_intern_2D Subroutine

  subroutine load_BC_intern_2D( tree, me, nSides )
    ! ---------------------------------------------------------------------------
    type(treelmesh_type), intent(in) :: tree
    !> Boundary condition construct to load the data into
    type(tem_BC_prop_type), intent(inout) :: me
    !> Required sides to set, defaults to 26.
    integer, intent(in), optional :: nSides
    ! ---------------------------------------------------------------------------
    integer :: iElem
    integer :: neigh_coord(4)
    integer :: my_coord(4)
    integer(kind=long_k) :: tID, neighID
    ! ---------------------------------------------------------------------------

    if (present(nSides)) then
      me%nSides = nSides
    else
      me%nSides = 26
    end if

    me%nBCtypes = 0

    allocate(me%BC_label(me%nBCtypes))
    allocate(me%boundary_ID(me%nSides, me%Property%nElems))
    me%boundary_ID = 0_long_k
    do iElem=1,me%Property%nElems
      tID = tree%treeID(me%Property%ElemID(iElem))
      me%boundary_ID(q__B, iElem) = -tID
      me%boundary_ID(q__T, iElem) = -tID
      if (me%nSides > 6) then
        my_coord = tem_coordofID(tID)
        ! Neighbor in the west
        neigh_coord = my_coord + [-1, 0, 0, 0]
        neighID = tem_IdOfCoord(neigh_coord)
        me%boundary_ID(q_BW, iElem) = -neighID
        me%boundary_ID(q_TW, iElem) = -neighID
        ! Neighbor in the east
        neigh_coord = my_coord + [+1, 0, 0, 0]
        neighID = tem_IdOfCoord(neigh_coord)
        me%boundary_ID(q_BE, iElem) = -neighID
        me%boundary_ID(q_TE, iElem) = -neighID
        ! Neighbor in the south
        neigh_coord = my_coord + [0, -1, 0, 0]
        neighID = tem_IdOfCoord(neigh_coord)
        me%boundary_ID(q_BS, iElem) = -neighID
        me%boundary_ID(q_TS, iElem) = -neighID
        ! Neighbor in the north
        neigh_coord = my_coord + [0, +1, 0, 0]
        neighID = tem_IdOfCoord(neigh_coord)
        me%boundary_ID(q_BN, iElem) = -neighID
        me%boundary_ID(q_TN, iElem) = -neighID
        if (me%nSides > 18) then
          ! Neighbor in the south-west
          neigh_coord = my_coord + [-1, -1, 0, 0]
          neighID = tem_IdOfCoord(neigh_coord)
          me%boundary_ID(qBSW, iElem) = -neighID
          me%boundary_ID(qTSW, iElem) = -neighID
          ! Neighbor in the south-east
          neigh_coord = my_coord + [+1, -1, 0, 0]
          neighID = tem_IdOfCoord(neigh_coord)
          me%boundary_ID(qBSE, iElem) = -neighID
          me%boundary_ID(qTSE, iElem) = -neighID
          ! Neighbor in the north-west
          neigh_coord = my_coord + [-1, +1, 0, 0]
          neighID = tem_IdOfCoord(neigh_coord)
          me%boundary_ID(qBNW, iElem) = -neighID
          me%boundary_ID(qTNW, iElem) = -neighID
          ! Neighbor in the north-east
          neigh_coord = my_coord + [+1, +1, 0, 0]
          neighID = tem_IdOfCoord(neigh_coord)
          me%boundary_ID(qBNE, iElem) = -neighID
          me%boundary_ID(qTNE, iElem) = -neighID
        end if
      end if
    end do

  end subroutine load_BC_intern_2D